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Abstract 


A linearized Euler solver for calculating unsteady flows in turbomachinery blade 
rows due to both incident gusts and blade motion is presented. The model accounts 
for blade loading, blade geometry, shock motion, and wake motion. Assuming that 
the unsteadiness in the flow is small relative to the nonlinear mean solution, the 
unsteady Euler equations can be linearized about the mean flow. This yields a set 
of linear variable coefficient equations that describe the small amplitude harmonic 
motion of the fluid. These linear equations are then discretized on a computational 
grid and solved using standard numerical techniques. For transonic flows, however, 
one must use a linear discretization which is a conservative linearization of the non- 
linear discretized Euler equations to ensure that shock impulse loads are accurately 
captured. Other important features of this analysis include a continuously deform- 
ing grid which eliminates extrapolation errors and hence, increases accuracy, and a 
new numerically exact, nonreflecting far-field boundary condition treatment based on 
an eigenanalysis of the discretized equations. Computational results are presented 
which demonstrate the computational accuracy and efficiency of the method and 
demonstrate the effectiveness of the deforming grid, far-field nonreflecting boundary 
conditions, and shock capturing techniques. A comparison of the present unsteady 
flow predictions to other numerical, semi-analytical, and experimental methods shows 
excellent agreement. In addition, the linearized Euler method presented requires one 
to two orders-of-magnitude less computational time than traditional time-marching 
techniques maJcing the present method a viable design tool for aeroelastic analyses. 
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Chapter 1 
Introduction 


In order to understand unsteady aerodynamic phenomena in turbomachinery, aeroe- 
IcLsticians require accurate and efficient models of the unsteady flow fields that result 
from blade motion and incident gusts. This is a formidable task since the flow is 
time-dependent, three-dimensional, compressible, and in general viscous. The objec- 
tive of this research is to develop a computational algorithm, based on the linearized 
Euler equations, that will accurately and efficiently predict unsteady subsonic and 
transonic flows in cascades. Specifically, the goal of the present research is to develop 
an algorithm which requires about the same amount of computer time to solve a sin- 
gle unsteady flow problem as is required to solve the equivalent steady flow problem, 
and which computes the unsteady flow field about as accurately as its steady flow 
counterpart. 

In the present study, particular emphasis is placed on the modeling of the resultant 
forces associated with unsteady shock motions encountered in transonic flows. It is 
important to model accurately these shock loads since the unsteady blade loading 
associated with the moving shock often accounts for about half of the total unsteady 
load. 

Another concern is the effective handling of the far-field computational bound- 
aries. Because the computational domain must be finite in extent, so-called far-field 
boundary conditions must be imposed on the far-field boundaries. If not properly for- 
mulated, however, these boundary conditions can reflect unsteady disturbances back 
into the computational domain, corrupting the solution. This can lead to incorrect 
predictions of the unsteady surface loading and hence incorrect conclusions as to the 
stability or response of the airfoil. In this report, a new exact far-field boundary con- 
dition formulation is presented. The method has the advantage that it is generic, and 
can be applied to a wide range of flow models (potential equation, Euler equations, 
Navier-Stokes equations) and can be extended to fully three dimensional analyses. 

1.1 Historical Development 

The techniques which have been developed over the last 30 years to analyze unsteady 
flows in cascades can be classified into one of the following categories: analytical and 
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semi-analytical methods or field methods. Using analytical approaches, the governing 
equations are solved exactly or through the use of point singularities. The complexity 
of the analytical problem requires simplifying assumptions which mahe the problem 
tractable. These assumptions, e.g., thin airfoils, no camber and light loading, severely 
limit the applicability of these methods since these assumptions are violated by most 
fan, compressor and turbine blades. 

One of the first classical analyses of unsteady flows in turbomachinery was de- 
veloped by Whitehead [52], who considered two-dimensional, incompressible flow 
through a cascade of flat plate airfoils. Since the steady flow was assumed to be 
uniform, there could be no accounting for the steady loading on the blade. White- 
head later extended his model to included the effects of steady loading and showed 
that steady loading plays an important role in bending flutter [53]. Later, Atassi 
and Akal [3] developed an inviscid incompressible flow model which included the ef- 
fects of camber, thickness, and steady blade loading. With this model they studied 
steady and unsteady incompressible flows through a cascade of thick, cambered air- 
foils. Their results, like Whitehead’s, showed the importance of steady blade loading 
on the unsteady flow in the cascade. 

The effects of compressibility were next investigated by Whitehead [54] and Smith 
[43]. With compressibility, the added complexity of acoustic modes and acoustic 
resonajice is introduced. Smith found that for unloaded cascades, the theory was 
successful in predicting the cut-off behavior as well as the amplitude of the acoustic 
waves. For steadily loaded cascades, however, the amplitudes of the acoustic waves 
were not predicted well. 

Several researchers [1,20,47] have investigated cascades of vibrating flat plates in 
supersonic flow which is axially subsonic. The results obtained using these models 
indicate that bending flutter will not occur at the reduced frequencies and Mach 
numbers at which real compressors actually exhibit bending flutter. Bendiksen [7] 
later used a perturbation analysis to include effects due to steady loading, thickness, 
camber and shock motion and demonstrated the important role of shock motion in 
flutter prediction. 

Although the assumptions invoked by these analytical and semi-analytical meth- 
ods limit analyses to low cambered and lightly loaded airfoils, significant insight was 
gained through the use of these elegant, albeit simplified, models. Specifically, to 
model accurately the flow in a caiscade, the effects of steady blade loading are critical. 
This means an effective model must contain the effects of camber and thickness, ais 
well as any discontinuities the flow field may permit, i.e. shocks. The evolution of 
computer capabilities, coupled with the cost and difficulty in experimentally evaluat- 
ing the unsteady aerodynamic performance of turbomachinery blade rows, h<is pro- 
vided an impetus for the development of increasingly more sophisticated field models 
which capture more of the underlying physics of the problem such as arbitrary blade 
geometries and complicated shock structures. These models, classified as field meth- 
ods, can be further separated into two groups, time marching and linearized methods 
(sometimes called time linearized). 

The most direct approach taken has been to develop time-accurate time-marching 
Euler [14,18,28,31] and Navier-Stokes [5,17,40,41] solvers capable of analyzing the 
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unsteady flows found in turbomachinery. This approach has merit in that the nu- 
merical implementation is relatively straightforward and complicated two- and three- 
dimensional nonlinear flows can be analyzed. There is a cost associated with this 
approach, however. Because of the large number of grid points and the requirement 
that the analyses be both accurate and stable, the maximum allowable time step used 
in the calculations must be small making these calculations prohibitively expensive 
for routine design use. 

In a linearized analysis, the flow is assumed to be composed of a nonlinear mean 
(or steady) flow plus an unsteady, harmonically varying small perturbation flow. 
Since the resulting unsteady flow field can be described by linear, variable coefficient, 
differential equations in which the explicit time dependence has been removed (that 
is the equations are solved in the frequency domain), the solution of the equations 
requires significantly less computational effort than that of the full nonlinear unsteady 
flow equations. Verdon et al. [47,48,49] and Whitehead and Grant [56] pioneered the 
development of linearized potential solvers as applied to the subsonic blade motion 
problem (the flutter or aerodynamic damping problem). Later, Hall and Verdon [26] 
and Caruthers and Dalton [11] extended this technique to handle incident vortical and 
entropic gusts using the flow decomposition technique of Atassi and Grzedzinski [4]. 
However, because of the inherent assumption that the steady flow be irrotational and 
isentropic, potential analyses are applicable only in the subsonic and low transonic 
flow regimes. In addition, three-dimensional flows with inlet swirl, typical of those 
found in turbomachinery applications, cannot be modeled using potential techniques. 

To overcome some of the limitations of the potential methods, investigators have 
begun to develop linearized Euler techniques. Ni and Sisto [39] proposed a technique 
whereby the linearized harmonic Euler equations were solved using traditional time- 
marching techniques by introducing a pseudo-time dependence. In this work, the 
flow was assumed to be isentropic and the results were limited to flat plate cascades. 
Hall and Crawley [24] later developed a direct method of solving the linearized Euler 
equations for realistic compressor designs, as well as for transonic channel flows. For 
transonic flows, they used shock fitting to model the motion of the shock. Although 
shock fitting does provide sharp, clearly defined shocks and shock motions, it is cum- 
bersome and complex and therefore less desirable than shock capturing. In a recent 
paper, Lindquist and Giles [35] outlined the conditions under which shock capturing 
correctly models shock motion in a linearized Euler analysis. The linearized Euler 
technique has shown the potential for dramatically reducing the computational cost 
required to solve unsteady subsonic [22,25] <ind transonic [23] flow problems in both 
two and three dimensions while still modeling the dominant physics of the unsteady 
flow problem. However, the method is still somewhat immature, and a number of 
modeling and computational issues require further research to make the method a 
viable technique. 
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1.2 Overview of Strategy Employed 

A linear cascade is an array of blades similar to a blade row in a turbomachine, 
except that it is usually two-dimensionaJ, as shown by Fig. 1.1. Calculations from 
two-dimensional cascades can be directly useful for rotating machines, provided that 
the hub- to-tip ratio is nearly unity, though corrections are often made for radial 
effects. The variables which define the cascade geometry are the blade shape, solidity 
(the ratio of chord-to-gap) and stagger angle (see Fig. 1.1). In the absence of unsteady 
excitation, the blades are identical in shape, equally spaced and their chord lines are 
oriented at the same stagger angle relative to the axial flow direction. 




Figure 1.1: Geometry of a rectilinear cascade 

Although viscous forces produce the wakes which are a prime source of forced 
response in turbomachinery, inviscid effects dominate the resulting wake interaction 
with a downstream blade row, subject to the limiting assumptions that the flow re- 
mains attached and that the boundary layers are thin. In the present analysis, there- 
fore, the fluid flow through a single cascade of airfoils is modeled as two-dimensional, 
adiabatic, and inviscid. Subject to these constraints, the governing equations of the 
fluid motion are the two-dimensional Euler equations. These equations are nonlinear 
and unsteady and thus generally difficult to solve efficiently. 

To simplify the solution process, the present method makes use of the fact that 
for small disturbances, such as are typically found in turbomachine applications, the 
unsteady perturbation flow is small compared with the mean flow and can reasonably 
be linearized about the mean flow. The resulting equations axe linear, variable coef- 
ficient, differential equations. By introducing a pseudo time dependence into these 
equations, the equations may be solved by time marching them to steady state. Since 
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the equations need not be marched time accurately, acceleration techniques such as 
multiple- grid acceleration and local time stepping may be employed. As previously 
indicated, along with efficiency, accuracy of the solution is important. 

The present analysis uses two procedures developed to improve the accuracy of 
the method. First is the implementation of a continuously deforming grid. While 
several nonlinear time-marching Euler codes have been developed that use moving 
or deforming computational grids [6,8,42,46], most previous linearized potential and 
Euler solvers have used computational grids fixed in space for solution of the lin- 
earized equations [24,48,49,50]. To calculate the flutter stability of an airfoil, one 
assumes the blades vibrate with a fixed amplitude, frequency, and interblade phase 
angle. The linearized equations are then solved to determine the unsteady pressure 
on the airfoil surface. However, since the blades move through the stationary grid, 
the boundary conditions that apply at the instantaneous location of the airfoil must 
be extrapolated back to the mean location of the airfoil where the boundary condi- 
tions are actually applied. Similarly, once the unsteady flow solution has been found 
on the fixed grid, terms involving the mean flow pressure gradient must be added to 
the unsteady pressure at the mean blade location to obtain the unsteady pressure at 
the instantaneous position of the blade. These gradient terms are difficult to evalu- 
ate accurately in practice, especially near the leading and trailing edges of fan and 
compressor blades. The resulting errors introduced at the airfoil boundary make the 
computed solution of the unsteady aerodynamic loads on moving airfoils inaccurate 
and sensitive to small errors in the computed steady flow field. 

To overcome some of the difficulties associated with fixed grid potential calcu- 
lations, Whitehead and Grant [56] introduced a linearized potential transformation 
that can be viewed as a rigid-body motion of the grid. Because the grid in effect 
moves with the airfoil, the scheme does not require extrapolation terms at the airfoil 
boundaries. While this approach allows rigid-body motions of airfoils to be analyzed, 
flexible mode shapes, which are common in turbomachinery aeroelasticity, cannot 
be analyzed. Recently, Hall [21] developed a linearized potential analysis that uses a 
deforming grid and is capable of analyzing both rigid-body and elastic airfoil motions. 

The present research presents a deforming grid linearized Euler solver that is 
suitable for the aerodynamic damping calculation of turbomachinery blades. It is 
demonstrated that the use of a deforming grid results in significantly more accurate 
predictions of the unsteady flow field and aerodynamic loads. Because the grid de- 
forms elcistically, the method can analyze both rigid-body and flexible blade motions. 
Furthermore, as implemented in the present method, the additional computational 
ancdysis required to use the deforming grid need only be calculated once before the 
start of the basic numerical integration scheme, hence the moving grid adds virtually 
no computational expense. 

The second gain in accuracy results from the implementation of numerically exact 
nonrefiecting far-field boundary conditions. Accurate unsteady flow computations 
demand nonreflecting boundaries to prevent any spurious reflections from corrupting 
the solution. Previous investigators [19,24,49,56] have derived the exact analytical 
behavior of linearized potential and Euler equations in the far-field and then matched 
these results to the numerical solutions computed in the near field. The present work 
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presents an alternative approach, whereby, the numerically exact, far-field behavior 
is found from an eigenanalysis of the discretized equations. The resulting eigenmodes 
are then used to construct the numerically exact, nonreflecting boundary conditions. 

1.3 Outline 

In Chapter 2, the governing nonlinear Euler equations are presented. The small 
disturbance assumption is introduced resulting in the linearized Euler equations which 
model the small disturbance unsteady flow field. Also introduced is a deforming 
computational grid which conforms to the motion of the vibrating airfoils. Finally, 
the analytical problem is closed with the introduction of the near-field boundary 
conditions. 

In Chapter 3, the basic (i.e., subsonic) numerical integration scheme is developed. 
The properties of the scheme, including accuracy, consistency and stability, are ex- 
amined using a model equation. Also discussed are grid generation techniques, and 
acceleration techniques such as multiple grid acceleration and local time stepping. 

In Chapter 4, a new numerically exact far-field boundary condition formulation 
is presented. To put this new technique into perspective, more traditional analytical 
boundary condition treatments are first reviewed. The new exact formulation, which 
makes use of the eigenmodes of the discretized small disturbance equations, is then 
presented. 

Chapter 5 presents numerical predictions of unsteady flows computed using the 
present method for several cascade configurations operating in the subsonic regime. 
Solutions obtained are compared to several approaches including analytical and semi- 
analytical methods, computational methods, and experimental data. Both the accu- 
racy and efficiency of the present method are examined. 

Chapter 6 introduces the theory and resulting modifications to the previously 
developed algorithm needed to model transonic flows using shock capturing. The 
idea of a conservative linearization is developed and the effects of nonconservation on 
the unsteady solution are explored. Test cases are then presented which depend on 
these modifications and the results are compared to nonlinear computational data. 

FinaJly, Chapter 7 concludes the report with observations about the usefulness 
of the present method as a design tool. In addition, conclusions from the present 
analysis and recommendations for future work in this area are presented. 



Chapter 2 

Governing Equations 


This chapter begins with the development of the integral form of the Euler equations 
valid for a deforming control volume. Next, the Euler equations are linearized to 
obtain a description of the small disturbance unsteady flows in cascades. First, a 
deforming computational grid which conforms to the motion of vibrating airfoils is 
defined. Then the Euler equations are linearized in a frame of reference attached to 
the moving grid. Subsequently, the preceding results are coupled together to obtain 
the linearized Euler equations valid on a continuously deforming grid. Finally, the 
analytic problem is completed by specifying appropriate near-field boundary condi- 
tions. 


2.1 The Nonlinear Euler Equations 

The general technique for obtaining the equations governing fluid motion is to consider 
a small, deformable control volume (Fig. 2.1) through which the fluid moves, and to 
require that mass, momentum and energy be conserved. This produces four equations 
which, when combined with an equation of state, provide suflUcient information for 
the determination of the four conservation variables: p, pit, pu, e, where p is the 
density, e is the internal energy, and u and v are the x- and ^-components of the fluid 
velocity. 

The present analysis makes the following assumptions. The flow is inviscid, im- 
plying that shear stresses are absent and hence the tangential velocity of the fluid at 
a solid surface boundary does not necessarily vanish as it must for a viscous fluid. A 
consequence of this simplification is that only surface forces which act perpendicular 
to the control surface (i.e., pressure) are present. In addition, it is assumed that there 
are no body forces present. Hence gravitational, electrostatic, ajid magnetic effects 
axe neglected. Finally, the flow is assumed to be adiabatic and the fluid is assumed 
to be thermally perfect. The consequence of the first assumption is that no heat 
is transferred to or from the surroundings. The second assumption implies the gas 
obeys the ideal gas law. With these simplifying assumptions, one can construct the 
governing equations of the fluid motion for a deforming control volume. 
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Figure 2.1: Typical control volume 

2.1.1 Continuity Equation 

The conservation of mass requires that the time rate of change of mass within the 
control volume be equal to the mass flux crossing the control surface. In integral 
form, 

i iL ■ £ r ■ ^^y^ = ° 

where / and g correspond to the displacement of the control volume in the x- and 
t/-directions. 


2.1.2 Momentum Equation 

Newton’s second law of motion states that the time rate of change of linear momentum 
of a given set of fluid particles is equal to the sum of the forces acting on the fluid. 
Using Reynolds’ transport equation, the conservation of linear momentum is given 
by 


^ jj pu dxdy + {^pv} + p - P‘^~^ dy - ^ (^puv - dx = 0 

^ II ^puv - ~ + P~ dx = 0 


where p is the static pressure. 


( 2 . 2 ) 

(2.3) 
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2.1.3 Energy Equation 

The first law of thermodynamics states (neglecting potential energy) that the time 
rate of change of internal energy plus kinetic energy of a system is equal to the rate of 
heat transfer into the system less the rate of work done by the system. Coupled with 
the assumption that the fluid is an ideal gas, and again using Reynolds’ transport 
equation, the conservation of energy is given by 


£ 

dt 


jj edxdy + {^puh^ - S^dy - 




dx = 0 


(2.4) 


where is the total enthalpy given by 


e + p 
P 


7 P 

7 - 1 p 


+ ^ (t/2 + v^) 


(2.5) 


2.1.4 Summary 

Using the results we have just obtained, we can now state in vector and conservation 
form the governing equations of the fluid. For a Cartesian coordinate system the 
Euler equations can be written as 

i IL ^ ^ £ (" - -Li^- 

where U is the vector of conservation variables and F and G represent the flux vectors 
described below. 
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It can be shown that the Euler equations in differential form are given by 

gU dF dG 
dt ^ dx ^ dy 


(2.7) 


(2.8) 


2.2 Linearization Of The Euler Equations 

Having developed the governing equations, in this section the nonlinear unsteady 
flow is decomposed into a nonlinear mean (steady) flow and a linear, harmonically 
varying perturbation flow. As previously stated, this analysis is performed on a con- 
tinuously deforming grid to improve the accuracy of the method. When the formal 
development is presented, the ensuing linearized Euler equations will contain inhomo- 
geneous terms which arise from the grid deformation. Hence, the system of equations 
to be discretized are linear, variable coefficient, inhomogeneous partial differential 
equations. 
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2.2.1 Computational Coordinate System 

Most previous linearized analyses have used computational grids which are fixed in 
space. For flutter calculations this implies that the cascade of airfoils vibrate through 
the stationary grid. To apply the appropriate upwash boundary condition at the air- 
foil surface, an extrapolation from the mean position of the airfoil to the instantaneous 
airfoil position must be performed. This procedure produces additional extrapolation 
terms in the airfoil upwash boundary condition; these extrapolation terms are quite 
large and difficult to evaluate accurately, especially in the regions near the leading 
and trailing edges. In this investigation, we avoid this difficulty by using a deformable 
grid which conforms to the motion of the airfoils. Since the airfoil motion is restricted 
to small amplitude, harmonic motion, the grid motion is similarly constrained. When 
implemented, no troublesome extrapolation terms are required in the application of 
the upwash boundary condition. 

The physical coordinates (x, y, t) are related to the computational coordinates (^, 
r;, r) by the following transformation. 


x((,ri,T) 


(2.9) 


= V + 

(2.10) 


= r 

(2.11) 


The functions / and g are chosen so that the motion of the grid on the boundaries 
conforms to the motion of the airfoils and are in general complex to include phase 
differences between adjacent fiirfoils. In the interior of the computational domain, 
it is desirable that / and g be continuous and smooth functions to minimize the 
truncation error of the numerical integration procedure. Thus, a natural choice is to 
let the values of / and g be determined by the solution of Laplace’s equation with 
appropriate boundary conditions. 

A typical example of unsteady grid motion can be seen in Fig. 2.2. These figures 
depict the computational and physical grids for the Tenth Standard Configuration 
airfoils undergoing a pitching motion with an interblade phase angle, o’, of 90®. In 
the computational coordinate system, the airfoils and the grid appear to be stationary; 
in the physical domain, the airfoils and grid continuously deform. 

2.2.2 Flow Decomposition 

Now that the grid motion has been defined, a similar expansion of the flow variables 
using the same assumption of small harmonic perturbations is proposed, i.e., 

U(4, g, t) = U(^, rj) -H u'(^, (2.12) 

Here, U is the vector of conservation variables representing the zeroth order or mean 
flow field, and u' is the vector of small perturbation harmonic amplitudes of the 
conser’vation variables representing the first order or unsteady flow field. It is now 
apparent that there are two sources on unsteadiness. These are: the unsteadiness 
associated with the small perturbation of conservation variables and the unsteadiness 



16 


CHAPTER 2. GOVERNING EQUATIONS 



Figure 2.2: Left: Unsteady grid in computational coordinate system t]). Right: 
Unsteady grid in physical coordinate system (z, y) for the case of a cascade of airfoils 
pitching about their midchords with an interblade phase angle, <t, of 90°. Flow 
calculations are performed using a single blade passage. Multiple passages are shown 
for clarity. 
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associated with the motion of the steady flow as seen by an observer attached to a 
frame of reference in the deforming computational grid. 

In a similar manner the flux vectors and the individual components of the state 
vector can be expanded in a perturbation series as follows: 


F(e,T 7 ,r) = F(U)+^uV‘-- (2.13) 

glQ 

= G(U)+^u'ei- (2.14) 


pu(tT},r) = A pu'(^,T])e^^^ 

pv((,r},T) = + pv'{(,T])e^‘^ 

e(l,7?,T) = e(^,?7) + e'(^,7/)e^‘^^ 


where dF/dV and dG/dV are the Jacobians, which can be written cis 
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(2.15) 


(2.16) 


(2.17) 


2.3 Linearized Euler Equations In Computational 
Domain 

Recall (Eq. 2.6) that the vector form of the governing equations for a deforming 
control volume is given by 

s IL ^ i ” 
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Substitution of the perturbation assumptions into the conservation equations de- 
scribed by Eq. 2.6 yields 


-f 


^ [dr] + dge^‘^‘^) 


-/ 

JC3 


/ 

^ \ 

G -f {d^ -I- dfe^'^) = 0 

aU / 


(2.18) 


Grouping terms of zeroth and first order results in the mean flow and small distur- 
bance equations respectively. The mean flow Euler equations represented in integral 
form are 

<bFdTj-<t Gd( = 0 (2.19) 

J C5 Jc9 

Equation 2.19 is the familiar integral representation of the two-dimensional Euler 
equations. 

Similarly, the resultant linearized Euler equations in integral form are 

ff ju>u'd^dT}+ ff jujVd^dg + ff juVdgdf 
JJcv JJcv JJcv 


- 1 - 


jf ^dg + ^u'dg-jujflJdg 


= 0 


( 2 . 20 ) 


Written in this manner, each first order term from Eq. 2.18 has an obvious corre- 
sponding term in Eq. 2.20. For example, the first row of Eq. 2.20 represents the first 
order perturbation of the first row of Eq. 2.18. In the first row of Eq. 2.20, the first 
term represents the unsteadiness in the first term of Eq. 2.18 due to the perturbation 
u* while the second and third terms in Eq. 2.20 represent the unsteadiness due to the 
grid motion. Rewriting Eq. 2.20 so that the homogeneous terms and inhomogeneous 
terms are separated gives 

ju, £ u' = 

-ju JJ U {d(dg -t- dfdr]) 

+ {fUdg - gVd() - i (Edg -'Gdf) 

Jc3 J CS 


( 2 , 21 ) 
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where the homogeneous terms (homogeneous in u') appear on the left and the inho- 
mogeneous terms appear on the right. 

It is apparent that the mean flow as described by Eq. 2.19 is independent of the 
unsteady perturbation flow. The perturbation flow, on the other hand, depends upon 
the steady flow in two ways. Recall that the linearized Euler equations are linear, 
inhomogeneous, variable coefficient differential equations. The variable coefficients 
(essentially the Jacobians) are functions of the mean flow field. Second, the inhomo- 
geneous portion of Eq. 2.21 (the right-hand side) is a function of only the steady flow 
field and the prescribed grid motion (for the gust response problem, the inhomoge- 
neous term is identically zero since the computational grid is non deforming). Hence, 
the right hand side of Eq. 2.21 need only be calculated once at the beginning of the 
iteration procedure. For this reeison, the use of a deforming grid adds very little to the 
total computational time of the algorithm. In addition, the Jacobians need only be 
calculated once and then can be stored for later use significantly reducing the number 
of required floating point operations per iteration. 

For numerical integration of Eq. 2.21, it is convenient to make the mean flow 
variables, U, and the perturbation variables, u', artificially time dependent. As 
suggested by Ni and Sisto [30], let 

U(^, T], r) = U(^, Tf, r) -f u'(^, 7], r)eJ"^ (2.22) 


Substitution of this assumption into the nonlinear Euler equations, Eq. 2.6, and 
again collecting terms of equal order gives 


_d 

dr 


// U d^drj + Fdr]- Gd^ = 0 

JJcv J cs Jcs 


(2.23) 


i IL li - !§"'''«) = 

-jcj ff \J {d^dg + dfdr]) + ju) i {fUdg - gVd() - <l (Fdg-Gdf) (2.24) 

JJ CV Jc3 Jcs 


Now both the nonlinear mean flow equations and the linearized Euler equations con- 
tain explicit time derivative terms making the equations hyperbolic in time. This 
allows solution of the equations using conventional time-marching algorithms. The 
equations are marched in time until the conservation variables reach their steady state 
values. Hence, the time derivative terms in Eqs. 2.23 and 2.24 are driven to zero and 
the original steady nonlinear and linearized unsteady Euler equations are recovered. 
Furthermore, since only the steady state values of the mean and perturbation flow 
conservation variables are desired, there is no need to time march the equations time 
accurately; multiple-grid and local time stepping acceleration techniques are used 
greatly reducing the computational time required to solve an unsteady flow problem. 

To obtain the corresponding differential form of the mean and linearized conserva- 
tion statements ( 2.23 and 2.24), we must relate the spatial and temporal derivatives 



20 


CHAPTER 2. GOVERNING EQUATIONS 


in the physical coordinate system to those in the computational coordinate system. 
This relationship is given by 



r ^ ^ ^ 

dx dx dx 




d_ 

dy 


dt dri dr 

dy dy dy 


dy 




I ^ 2 ji Sz 

L dt dt dt 




(2.25) 


where [J] is the Jacobian matrix associated with the coordinate transformation. A 
similar transformation matrix can be obtained for the inverse problem, and will be 
denoted [J]~^ as shown below 



dx dy dt 
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dx dy dt 

dy dy dy 


dy 


[J]-^ 


dy 



dx dy dt 

L dr dr dr J 



VI/ 


(2.26) 


Substituting the small motion assumption of Eqs. 2.9 through 2.11 into the [J]“^ 
matrix and inverting, the [J] matrix is found to be given by 


1 - /« -9i 0 


[J] = 


-fy ^-9y 0 

-fr -9r 1 


(2.27) 


where only zeroth and first order terms have been retained. 

The derivatives in the physical coordinates can now be represented in terms of 
the grid coordinates as follows: 


d _ d ( d d\ 
dx~ di 

dy dy \ "' d^ ^ dy ) 

I.'] 

dt~dr y^d^'^^^dy) 


(2.28) 

(2.29) 

(2.30) 


Thus, in vector form, the differential representation of the linearized Euler equations 
in the computational domain is given by 


juu + — + 


^ ^.5u' - 




d^ dV dy dV 
d^ ^ dy) didy) \dy d^ dydy) 



2.4. NEAR-FIELD BOUNDARY CONDITIONS 


21 



Figure 2.3: Clcissification of boundary types in computational domain 

2.4 Near-Field Boundary Conditions 

To complete the specification of the unsteady flow problem, boundary conditions must 
be specified around the entire computational domain as shown in Fig. 2.3. There are 
essentially three types of unsteady boundary conditions. The solid surface, or no 
through flow condition dictates that there be no flux of mass through solid surfaces. 
Periodic boundary conditions are applied upstream and downstream of the cascade 
to reduce the computational domain to a single blade passage. Finally, far-field 
boundary conditions are necessary to prevent spurious reflections of outgoing waves 
at the upstream and downstream computational boundaries (the development of the 
far-field boundary conditions is deferred until Chapter 4). 

2.4.1 Solid Surface 

The solid surface boundary condition is necessary to ensure that no flow passes 
through the surface of the airfoil. The development of this boundary condition will 
begin with the full unsteady nonlinezir representation of this statement which is given 
by 



22 


CHAPTER 2. GOVERNING EQUATIONS 



Figure 2.4: Illustration of position vectors for flow tangency boundary condition 


(2.32) 

Figure 2.4 illustrates the location of the airfoil surface at its mean and perturbed 
location. R(s, r) describes the location of the reference airfoil’s surface at time r, 
where s is the distance along the airfoil surface. 

In a manner analogous to that developed in the section pertaining to flow decom- 
position, the unit normal vector and the airfoil’s position vector can be expanded in 
perturbation series: 


R(s, r) = R(s) + r'(s)e-'‘*'’’ 

(2.33) 

n(s,r) = n(s) •+• 

(2.34) 

V(s,r) = V(s)-|-v'(s)e><^ 

(2.35) 

Substitution of Eqs. 2.33 through 2.35 into Eq. 2.32 yields 


(V + v')-(n + n')-^^^i^-(n + n') 

(2.36) 

where 


r' = {f,g) 

(2.37) 
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Recall that / and g are the functions that represent the first order unsteady grid 
motion. Collecting terms of zeroth and first order describes the steady and unsteady 
flow tangency conditions respectively, 

V-n = 0 (2.38) 

v'-n = —Y-n' + ju{f,g)-n (2.39) 

Equation 2.38 is clearly a statement restricting the steady flow to be tangent to 
the blade surface. The unsteady counterpart to the flow tangency condition, Eq. 2.39, 
states that the normal perturbation flow velocity must be equal to the upwash on the 
surface of the airfoil attributed to its unsteady motion. The first term on the right 
side of Eq. 2.39 is the upwash associated with a rotation of the blade surface, whereas 
the second term is the upwash associated with a translation of the blade surface. 
The benefit of using a deforming grid is now clearly evident. There are no terms in 
Eq. 2.39 that represent extrapolations of the blade’s position to a mean location in 
order to apply the unsteady boundary condition. Finally, it should be noted that the 
upwash is only a function of the steady flow and the specified grid motion. Both of 
these quantities are known prior to the iteration procedure. 

2.4.2 Periodicity 

Because the unsteady flow is governed by linear equations, the response of the 
cascade to blade motion or a gust may be found by decomposing the disturbance into 
a sum of travelling wave modes, each having a different frequency, a;, and interblade 
phase angle, <r. By superposition, the total response of the cascade is equivalent to 
the sum of the response to each of the individual modes. Therefore, without any 
loss in generality, only one travelling wave mode need be considered at a time. The 
unsteady periodicity condition upstream and downstream of the cascade is then 

V + G)= u'(^, r/)eJ^ (2.40) 

where G is the blade- to*blade gap (see Fig. 2.5). This periodicity condition allows 
the computational domain to be reduced to a single blade passage which significantly 
reduces memory requirements and computational time. 

2.5 Summary 

In this chapter, the linearized Euler equations were developed. Included in this anal- 
ysis was the effect of a deforming computational grid. Near-field boundary conditions 
were presented to complete the specification of the unsteady aerodynamic model. 
In the next chapter, a numerical integration scheme for solving the linearized Euler 
equations will be presented. 
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Figure 2.5: Illustration of the periodic boundary condition 



Chapter 3 

Basic Numerical Integration Scheme 


Over the last two decades numerous computational methods have been developed 
to solve the steady Euler equations (e.g., [13,37,39]). One goal of the present re- 
search, is to adapt proven steady flow solvers to solve the unsteady linearized Euler 
equations. An advantage of this approach is the same algorithm can be used in the 
solution of both the steady, nonlinear Euler equations and the linearized, unsteady 
Euler equations. Additionally, the Lax-WendrofF algorithm method selected has been 
proven to be an efficient and effective way of solving the nonlinear Euler equations 
and is believed to be equally effective when applied to the linearized Euler equations. 
With this in mind, the main goal of this chapter is to present the basic numerical 
integration method used in the present analysis. In addition, issues of grid generation 
and numerical properties of the integration scheme such as stability, consistency, and 
accuracy are investigated. 

3.1 Unsteady Grid Generation 

The primary focus of the present work is the development of a numerical tool aeroelas- 
ticians can use to calculate unsteady flows in turbomachinery applications. Toward 
that end, the steady Euler and linearized Euler equations will be discretized and 
solved using a finite volume technique. Before proceeding, however, a computational 
mesh must be constructed to specify the discrete locations where the solution is to be 
found. It is well known that the quality of the grid used can significantly affect the re- 
sults. Therefore, while grid generation is not the primary focus of this research, a brief 
development of the grid generation techniques used is presented for completeness. 

The problem of grid generation has been an area of active research for many years 
[2,44]. Generally, the greater the resolution the more accurate the results. On the 
other hand, the greater the resolution, the more computational time and resources are 
required to solve the problem. This introduces one of the fundamental problems of 
grid generation: how to balance the resolution required to obtain a reasonable solution 
with an acceptable amount of computational resources. Another desire is that the 
grid lines be nearly orthogonal everywhere as most finite difference cind finite volume 
schemes produce large truncation errors when the grid is highly skewed. Although this 
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may be possible, it comes at the expense of added complexity in the grid generation 
algorithm. In addition, some skewedness is permissible without significantly affecting 
the solution. Finally, adjacent cells should be approximately the same size, since 
rapid changes in the size of the computational cell can introduce error as well as lead 
to stability problems. These factors and others axe taken into consideration when 
ascertaining the quality of the mesh. It is for these reasons that there is an art 
associated with grid generation that is outside the scope of the present research. 

Since the unsteady grid motion is comprised of a steady (mean) grid plus a har- 
monically varying small perturbation, it is necessary to compute first a steady grid, 
about which the unsteady grid motion can be linearized. In the generation of the 
steady grid, a modified version of an elliptic mesh generator developed by Thompson 
[44] is used. The principle behind this grid generation technique is to prescribe the 
boundary points in the physical dom«dn and then map the specified computational 
grid into the irregular physical domain by solving an elliptic PDF with the appropriate 
boundary conditions. 

With the steady grid generated it is now necessary to generate the unsteady grid 
motion. Recall that the position (x,y) of a grid point is given by 

a; = ^ + /(^, 

1/ = 7/ (3.1) 

For cases involving blade motion, the complex amplitude of the grid motion 
must be specified at every grid point in the computational domain. Since the mode 
shape of the blade vibration is known, the quantities / and g are prescribed on 
the solid surface boundaries to match the prescribed airfoil motion. Upstream and 
downstream of the airfoil, the values of / and g on the boundary are prescribed to 
“smoothly” vanish as the fax field is approached. The requirement that the motion 
of the grid vanish in the far field is made to simplify the far-field boundary condition 
development in Section 4.3. Additionally, the grid motion must satisfy a complex 
periodicity condition similar to that prescribed on the unsteady flow variables. 

With the values of / and g now prescribed around the entire computational do- 
mciin, the problem now is to determine the motion in the interior. In principle the 
distribution could be arbitrary. However, a smooth grid distribution reduces trun- 
cation errors in the unsteady flow computation. Therefore, in the present work, the 
motion of the grid on the interior is described by Laplace’s equation so that 


vv = o 


V^g = 0 (3.2) 

subject to the Dirichlet boundary conditions specified above. A finite element scheme 
is used to discretize the equations on the steady grid. The equations are then solved 
directly using LU decomposition. Because the motion of the grid is harmonic, the 
solution of the grid motion need only be determined once before the iterative solution 
procedure begins. Figure 2.2 shows the computational grid (steady) and physical grid 
(unsteady deforming) for a typical unsteady flow calculation. 
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3.2 Basic Lax-Wendroff Integration Scheme (One- 
Dimensional Model) 

When deciding which integration scheme to use in the present research, the following 
factors were considered. First, it is desirable to take advantage of previous experience 
so as not to duplicate elfort. Second, the same algorithm should be used in both the 
steady and unsteady codes. Third, the scheme should be computationally efficient 
making use of modern convergence acceleration techniques. Hence, a logical choice 
for the present research was an adaptation of the Ni scheme [13,37] which is itself a 
variant of the familiar Lax-Wendroff scheme. The Ni scheme, as traditionally applied 
to the steady nonlinear Euler equations, is marched in time until a converged solution 
is obtained. However, since the linearized Euler equations are cast in the frequency 
domain, they are not explicitly dependant on time. Hence, as suggested by Ni and 
Sisto, a pseudo time term must be introduced as discussed in section 2.3. The resulting 
linearized Euler equations in differential form are 


du' 

dr 


-f juu' + 


dU dri au 


(3.3) 


where b is the right side of Eq. 2.31. Equation 3.3 resembles Eq. 2.31 with the 
exception that an additional time dependent term now appears on the left side of the 
equation. This pseudo-time dependent term is driven to zero by time marching the 
modified linearized Euler equation to steady state, thereby recovering the solution to 
the original linearized Euler equations (Eq. 2.31). 

Ni’s Lax-Wendroff scheme can easily be applied to the linearized Euler equations 
with some minor modifications. To simplify the description of the method, the inte- 
gration scheme will be presented using the following one- dimensional model equation. 


du dF{u) 
dt dx 


= 0 


If Eq. 3.4 is linearized, the following expression is obtained. 


(3.4) 


. , d 



= 0 


(3.5) 


Introducing the pseudo time dependence results in the following expression for the 
perturbation equation. 


du' d 



= 0 


(3.6) 


This new model equation resembles the linearized Euler equations as written in Eq. 3.4 
and it is for this equation that a variation of Ni’s scheme will be presented. The 
variation from the original development arises from the introduction of the joju' term 
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Figure 3.1: Typical computational cell for one-dimensional Ni scheme 


which is a result of the linearization process. (In the following development the primes 
will be omitted for clarity) 

To begin, the solution u is expanded in a Taylor series, resulting in 


ti"+^ 


u" -1- Sui 



(3.7) 

(3.8) 


where the superscript n refers to the computationcd time level and the subscript i 
refers to the grid location. In Eq. 3.8 is the correction in the solution at the node 
from one time level to the next. Rewriting Eq. 3.6 gives 


du r. d fdF \ 


(3.9) 


Differentiating Eq. 3.9 with respect to time and applying the chain rule yields 



. du d 

(dFduX 


L 

ydUdt J 


(3.10) 


Finally, substituting Eqs. 3.9 and 3.10 into Eq. 3.8 and evaluating at grid index, z, 
and time level, n, with centered finite difference expressions for the spatial derivatives 
results in the following expression for the change in the state variable from one time 
level to the next. 


Sui = 


^ Au.-i + u.) (u;-bu,>i) 

2 V 2 2 


2Ax 


+ 


f: - F.' , 
2Ai 


At 
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t+l/2 



where 


■ du 


Ui 


(3.U) 


(3.12) 


One advantage of Ni’s scheme is that it is cell based. To demonstrate, we introduce 
the following cell based notation: 


aub = (f;:, - F,'”) 


At 

Ax 



An 


B 


i-l/2 



(3.13) 

(3.14) 


+1/2 


and 


^ (“ 7 - 1 + “r) 

= \ (“• + «7+i) 


(3.15) 


Here Aug and Au^; are called the changes at cells B and C respectively. The subscript 
B refers to the control volume lying between grid points i — l and i and, the subscript 
C refers to the control volume lying between points i and i + 1 as shown in Fig. 3.1. 
Cast in this manner, this implies the correction to point i is composed of two parts. 
The first part. 


««,)b = 5 


, At 

AUg + ~ 


jwAt o;2At2_ _ ( dF 

+ — — “B - “B 1 ^ 


At 

Ax 


(3.16) 


is the correction due to the change, Au^, taking place in control volume B. The 
second part, 


(«“i)c = I 


, At 

Auq — AF«— JuAuq 

^ Ax 


jwAt 


Auc + 


u;2At2 


■uc + jwAt uc ^ 



(3.17) 


is the correction associated with the change, Au,;;, occurring in control volume C. 
Note that the expressions for (^uJb and (^u,)^; are similar with the exception of the 
sign on two terms. Although these expressions resemble those presented by Ni [37], 
there are some important differences. Most notably, there are terms which have a 
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frequency dependence that are the result of the linearization process. These terms 
will be shown in a subsequent section to play a major role in the stability of the 
system. 

To summarize, a control volume centered integration method is devised that em- 
ploys the following operations in the determination of the correction to an arbitrary 
grid point. 

(1) Sweep through all cells computing the change in each cell using 


Auc + ju;uc = ^ (f;^ - + ju 


Ui + U.-+1 ' 


At 


(3.18) 


(2) Determine the corrections to point i and z-|-l through the use of the distribution 
formulae 




Aucr — ^F„- jojAuc 

^ Ax 


juAt _ 

- - Auc H — uc + ;wAt uc 


dF 

dU 


At 

Ax 


and 




, At 

Auq -1- ~ 


jujAt 


2 " ■ 2 
(3) Update the dependent variable by 


^ o;2At2_ . _ dF 

+ — 7^ — - JwAt Uc 


dU 


At 

Ax 


= u" -|- 6u{ 

t I * 

Sui = 

where both and (^u,)c are known from step 2. 

(4) Apply appropriate boundary conditions at the inlet and exit. 

(5) Continue the iteration procedure until the solution converges. 


(3.19) 


(3.20) 


(3.21) 

(3.22) 


3.3 Basic Lax-Wendroff Integration Scheme (Two- 
Dimensional Euler Equations) 

The bcisic integration procedure outlined for the one dimensional model equation in 
the previous section can be extended and applied to the linearized two-dimensional 
Euler equations on a non-orthogonal, curvilinear grid (see Fig. 3.2 for nomenclature). 

As in the development of the integration scheme for the one dimensional model 
equation, the following five steps are performed during each iteration of the numerical 
algorithm. For example, at each control volume, say C, the process is: 
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Figure 3.2; Nomenclature used in the description of the computational cells 


(1) Sweep through all cells computing the change in each cell using 


Auc + juAt uc 


AV 1 2 Vvnw 'Isw/ 2 VSntx; SswJ 


f +f' G' +G' 1 

“(in, - 1 „) - - ” ({n, - e„) 

f' +F' G' +G' 

^^ 1 .. - 1 ») - “ ( 4 . - 


+ \ (in. - Inn,) - ^ ( 4 „ ~ U)] } ( 3 - 23 ) 


— 2 ^nw)i7]ne Vsw') (Cne ^iui)(^je ^ntu)] 


(3.24) 


(2) Determine the corrections to the four corner nodes through the use of the 
distribution formulae 

1 juAt 

= 4 - ^fc - ^9c 

1 r jujAt 

(^«nw)c = j - A/c + Age 

1 juiAt 

(^^nc)c = 1 + ^fc + ^9C 
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{S^se)c = J 


juAt 

Auc + Afc - Age 


(3.25) 


Again, the following notation has been introduced to simplify the expressions 


A/c = ^ - AG'^A?') 

Asc = ^ (ACcAr - AF^A,”) (3.26) 

where ^ 

+ ^ne - Lw - 

2 Vae) 

A-f"* = - iinc + Le - Lw - Lw) 

Ag^ = i + Vac - Vaw ~ Vnw) (3.27) 


(3) Update the dependent variable by 

= u" + 8U; 

t I ' 

Su,- = (Sui)^ + (Su^)g + (Su,)e + (Sui)^ (3.28) 


where (5u,)^, (Su^Jg, (Su^)c, and (Sui)g known from step 2. 

(4) Apply appropriate boundary conditions at the inlet and exit. 

(5) Continue iteration procedure until the solution converges. 

It should be noted that instead of a single equation (as in the model problem) 
the Euler equations represent four separate equations. As such, the changes and 
corrections described in the previous development are now vector quantities described 
below. First, the changes in conservation variables are 


Au' = 


(Ap') 

(Apu') 

(Apv') 

(Ae') 


Now, AF' and AG' can be obtained by the following expressions, 


and 



(Apu') 

U (Apu') + U (pAu') + Af/ 
V {Apu') + U {pAv') 
h.{Apu') + U (pAh'J 



{Apu') 

U {Apv') + V (pAu') 

V {Apv') + V {pAv') + Ap' 
K {Apv') + V {pAh'J 


(3.29) 


(3.30) 


(3.31) 
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where 


{pAu') 

ipAv') 

{Aj/) 

{p^K) 


{Apu') - U(Ap') 
{Apv’) - V {Ap‘) 


1 


( AeO - - (£/ {Apu') + y {Apv') + U {pAu') + V {pAv')) 


(7-1) 

{Ae') + {Ap')-\{Aj/) 


(3.32) 


The final consideration is the time step, At. For speed in convergence, the larger 
the time step the better, however, the time step must not violate the CFL (Courant, 
Friedrichs, and Lewy) condition. For the scheme to be stable, the numerical scheme 
cannot propagate the solution faster than the physical wave speeds. This is the so 
called CFL condition. To this end, the following criteria is used to determine the 
maximum allowable time step. 


Ay 

\UAtj^ -VA(^\+ a Al 

where 

A/= + , 

and a is the speed of sound. 


Ay \ 

(3.33) 

\UAri”' — yA^'^l + a Am J 

Am = y^(A^”‘)^ + (Arj^)^ 

(3.34) 


At < min 


3.4 Smoothing 

The numerical equations presented can be applied to obtain solutions to inviscid, 
rotational (and/or irrotational), subsonic, transonic and supersonic flow problems. 
However, for subsonic flows a small amount of artificial viscosity must be added to 
suppress the spurious sawtooth solutions permitted by the numerical solution tech- 
nique. Furthermore, for transonic and supersonic applications the addition of viscous 
like terms are needed to stabilize and capture shocks. Although it has been deter- 
mined that some level of smoothing is imperative, care must be taken so as not to 
add so much as to diminish the accuracy of the numerical algorithm. In the present 
analysis, two different smoothing terms have been introduced. The first, used in trcin- 
sonic and supersonic flow regimes, is termed second difference smoothing. Although 
this smoothing is first order accurate it only need be applied in the regions near the 
shock. The second, used in both subsonic and transonic regimes, is termed /ourt/i dif- 
ference smoothing and is effective at reducing the numerical oscillations introduced. 
As implemented, the solution procedure with only fourth difference smoothing re- 
mains second order accurate. 

To illustrate how the smoothing is added to the corrections (previously deter- 
mined), we again consider the one-dimensional model problem. The correction to the 
ith node with smoothing is found to be 


Sui = {Sui)g + {Sui)c + cTj (ui+i - 2ui -f- u,_i) 
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+ <^4 (“i-2 - 4ui-i + 6u; - + W.+2) (3.35) 

where and cr^ are coefficients of smoothing indicating the level of second differ- 
ence and fourth difference smoothing added to the correction. Typically the values 
range from 0.0 to 0.05 for the fourth difference smoothing coefficient, and 0.0 to 
0.1 for the second difference smoothing coefficient, <T 2 - Equation 3.35 would seem 
to indicate that the fourth difference smoothing has a difference stencil of five grid 
points (corresponding to 25 grid points for two-dimensional Euler). As implemented, 
however, the fourth difference operator can be generated by taking the second dif- 
ference of the second difference operator. As a result, the difference stencil for the 
one-dimensional model equation remains 3 grid points (corresponding to 9 grid points 
for two-dimensional Euler) and hence the smoothing operators can be adapted to the 
cell based distribution formulas previously developed. 

3.5 Properties of Ni’s Scheme 

At this point in the analysis, the analytic problem developed in Chapter 2 has been 
discretized. The dependent variables are now defined at discrete locations or nodes. 
Spatial derivatives have been approximated using finite volume techniques resulting 
in a system of algebraic equations. Thus, the original problem involving a continuous 
domain and partial differential equations hcis been transformed into a discrete system 
of algebraic equations. At this time it is important to determine whether the solution 
obtained by solving the algebraic system is a good approximation to the original 
system of PDEs. Among the issues considered in the following sections are consistency 
and stability, 

3.5.1 Accuracy and Consistency 

Simply stated, the truncation error of a scheme can be defined as the difference 
between the partial differential equation and the finite difference approximation to it. 
In other words, 

ErTort^^^^^tion = PDE - FDE (3.36) 

The order of the scheme is a meaisure of the magnitude of the truncation error. 
Consistency deals with the extent to which the finite-difference equations approximate 
the partial differential equations. A finite difference representation of a PDE is said 
to be consistent if it can be shown that the truncation error vanishes as the mesh is 
refined. 

A formal development of these two concepts will now be presented. Again, it 
is convenient to illustrate these ideas with a less complicated model equation which 
resembles the linearized Euler equations in form and behavior. The model equation 
considered is the one- dimensional wave equation, which when linearized and reintro- 
ducing a pseudo time dependence takes the following form 


Uf + juu + CUj. = 0 


(3.37) 
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where the subscripts t and x represent differentiation with respect to time and the 
x-direction respectively. Note that Eq. 3.37 is identical to Eq. 3.6 when c is equivalent 
to the Jacobian {dFfdU). With this in mind, the finite difference expression for this 
model equation is then 


u""*"' = u'y — jot 


u'} , + 2u^ + u?, 

t— 1 » tH 


^) + 5 [y« - 1) - K-,) 


+ — { - 2u7 + , ) 

2 V ‘+1 ‘ ‘-V 2 V 


u^_, + 2u^ + ’ 


(3.38) 


where the following simplifying notation heis been introduced 


a = wAf 

A = c— 
Ax 


(3.39) 


where a is a reduced frequency based on the time step, and A is the so-called CFL 
number. Equation 3.38 is the final finite difference expression for the model equation. 
Using this equation and Taylor expanding about u^, (for clarity the subscripts i and 
n are dropped) the following modified equation is obtained. 


Uj -f juu + cu^ = u 
w^A^Ax^ 


.w^A^Ax^ w'^A^Ax^ 


-J 


+ 

+ 


2c 


+ i- 


6c^ 
cj^A^Ax^ 


24c3 


6c^ 


d" ^rr; 
“f" ^xxix 


Ax2A u;2 AAx3\1 

rj'^T ) 

cAx2 .o;Ax3 


6 

cAx3 

24 


(A-A3)-i 


6 

wAx^ 
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0)2 A Ax® 
96c 


(3.40) 


The modified equation is the partial differential equation which is actually solved 
when a finite difference method is applied to a PDE. It is important to emphasize 
that the equation obtained after substitution of the Taylor-series expansion must be 
used to eliminate the higher-order time derivatives rather than the original PDE. This 
is due to the fact that a solution of the original PDE does not in general satisfy the 
difference equation. 

Upon examination of Eq. 3.40 it is obvious that the left hand side of the equation 
is the original model PDE. The right hand side of the equation is known as the 
truncation error. It is important to note that as the grid is refined the error terms 
on the right hand side go to zero like Ax2. As implemented, the scheme is therefore 
second order accurate in both space and time and hence is consistent. Although this 
alone is not surprising from previous cinalyses on the Lax-Wendroff scheme, what is 
unique about this analysis is that the inclusion of the complex source term (jwu) 
which changes the character of the truncation error. Any term with an uj dependence 
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would traditionally not appear in a time domain analysis. As will be shown, these 
terms not only a,fFect the character of the solution, but have an important influence 
on the stability of the system. 


3.5.2 Stability 

Lax’s equivalence theorem states that given a properly posed initial value problem 
and a consistent finite-difference approximation to it, stability is the necessary condi- 
tion for convergence [2]. Convergence means that the solution to the finite- difference 
equation approaches the true solution to the PDE having the same initial and bound- 
ary conditions as the grid is refined. A stable numerical scheme is one for which errors 
from any source (round-off, truncation, etc.) do not grow as the calculation proceeds 
from one time step to the next. 

To illustrate these ideas, we consider again the difference expression for the one- 
dimensional model equation (the frequency domain version of the wave equation) 
given by Eq. 3.38. Two separate stability analyses are performed. The first analysis 
employs a classic Von-Neumann stability analysis. In this approach the domain is 
assumed to be infinite, hence no special consideration of the boundaries is taken into 
account. It will be shown that the modified Ni’s scheme, as applied to the model 
equation, is unconditionally unstable as determined for an infinite computational 
domain. Next, an analysis is presented which determines the stability of the modified 
Ni scheme on a finite domain. It will be shown that the finiteness of the computational 
domain has a stabilizing effect on the system and thus allows for converged solutions 
to the finite difference equation. 


Von-Neumann Stability Analysis 


Using the previously defined notation, Eq. 3.38 can be written as. 


= 


d- 

+ 




a 


2 V2 


1 


A 


( o + + o (1 + ” 




a 


a‘ 


1-J--A2 

•^2 4 


.ct 


“«+i 


-^2 U '^ 




2(l + ^)-y 


(3.41) 


To perform a Von-Neumann stability analysis, the solution of Eq. 3.41 is assumed to 
take the form 

u" = (3.42) 


where 


m sAx 
i 

»— 1 

yTi _ |jngjfcm(t+l)Ax 

>+l 


(3.43) 
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Substitution of Eqs. 3.42 and 3.43 into Eq. 3.41 gives 




= G = 




+ 

+ 


Q-jkmAx 




+ ^(1 + ^) 


1 



g+jfcmAx 




+ '^) 




(3.44) 


For stability, the magnitude of the amplification factor, 1G| , must be less than or 
equal to one for all wave numbers, fc„,, at the desired frequency, uf. Based on this 
criteria it would indicate that some restrictions might be placed on At and Ax or 
more generally the Courant number, A, (also known as the CFL number) for the 
system to be stable. 


A = c— 
Ax 


(3.45) 


Figure 3.3 shows that the stability limits of the one-dimensional wave equation in 
the time domain (a? = 0). As expected, the results show that CFL numbers greater 
than unity cause instabilities to occur. This agrees with classical analyses of the time 
domain wave equation [2]. Figure 3.4 demonstrates the effect the ju;u source term 
has on the system’s stability. In this example, the reduced frequency, k = a/A, is 
1.0. It is evident that |G| is greater than one for some wave numbers regardless of 
the Courant number indicating unconditional instability. However, the influence of 
the boundary conditions has not been included in this analysis. In a Von-Neumann 
stability analysis, the implicit assumption is that the computational domain is infinite 
(or that the computational domain is periodic). In the next section the effects of the 
boundaries on the system’s stability will be investigated. 
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Figure 3.4: Von-Neumann stability analysis, one-dimensional model equation: k 
1.0, A = 0.6, 0.8, 1.0, 1.2. 






3.5. PROPERTIES OF NFS SCHEME 


39 


Eigenvalue Stability Analysis 

In the previous section, a traditional Von-Neumann stability analysis indicated that 
the discretized version of the model equation was unconditionally unstable. Strictly 
speaking, however, the Von-Neumann analysis holds only for infinite or periodic com- 
putational domains. Furthermore, it is the longer wavelength modes that are most 
unstable. These are the modes most likely to be influenced by the finiteness of the 
computational domain. 

The finite difference representation of the model PDE is a so-called explicit scheme. 
An explicit scheme is one for which only one unknown appears in the difference 
equation in a manner which permits evaluation in terms of known quantities. In 
matrix notation the finite difference expression can be written as 

{u-+^} = [A]{u-} (3.46) 

where [A] has the following structure 


[A] = 


B.C. 

XXX 

XXX 

B.C. 


(3.47) 


It is apparent that the stability is now governed by the spectral radius of [A]. If the 
maximum eigenvalue of [A] is less than or equal to unity, the scheme is stable. 

Figure 3.5 shows the results of an eigenanalysis of the finite difference represen- 
tation of the one-dimensional, time domain (w = 0) model equation for various CFL 
numbers. For this figure, periodic boundary conditions were implemented in an at- 
tempt to duplicate the results determined by the Von-Neumann ainalysis shown in 
Fig 3.3. Note the excellent agreement between the Von-Neumann analysis and the 
eigenvalue analysis. Next, Fig. 3.6 provides a similar analysis for the unsteady case 
when the reduced frequency, fc, is 1.0. Again periodic boundary conditions were 
used and the results agree with those presented in Fig. 3.4 determined from the 
Von-Neumann analysis. These results indicate that the discretized unsteady repre- 
sentation of the one-dimensional wave equation are unconditionally unstable when a 
periodic computational domain is used. 

In order to mimic the behavior of the linearized Euler equations more closely, 
we now replace the periodic boundary conditions previously used with nonreflecting 
boundary conditions at the inflow and exit boundaries (these boundary conditions 
are analogous to those used in the linearized Euler analysis which will be presented in 
Chapter 4). Figure 3.7 shows the effect the boundary conditions have on the system’s 
stability. In this case, the eigenanalysis of the discretized system indicates that the 
Lax-Wendroff scheme is stable. In effect, the finiteness of the computational domain 
has a stabilizing effect on the numerical algorithm. Figure 3.7 presents the eigenvalues 
of the system for a range of reduced frequencies, k, for a CFL number. A, of 0.8. This 
was accomplished by varying the number of grid points (thereby changing Ax) in a 
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Figure 3.5: Eigenvalue analysis of one-dimensional model equation: k = 0.0, A = 0.6, 
0.8, 1.0, 1.2 (periodic boundaries) 

duct of length 1.0. The frequency based on the duct length was held constant, uj = 
5.0. The interesting behavior to note here is that the “root locus” determined by 
the eigenanalysis approaches that of the Von-Neumann analysis as the grid is refined. 
This is due to the fact that in the limit the domain appears to become infinite (Ax 
approaches zero) and hence more closely models the Von-Neumann analysis. 

The numerical algorithm used in the present research has been shown to be both 
consistent and stable as long as the boundary conditions used in the far field stabilize 
the long wavelength modes. By virtue of Lax’s equivalence theorem the necessary 
conditions for convergence have been met. 

3.5.3 Conservation 

The PDE’s of interest in this research all have their basis in physical laws such as the 
conservation of mass, momentum and energy. Such a PDE represents a conservation 
statement at a particular point. The criteria necessary for a finite difference scheme 
to be considered conservative is now discussed. Consider the integral form of the 
continuity equation for a fixed control volume 

f V • pYdV = <f pY-ndS = 0 (3.48) 

Jcv Jc3 

To determine whether the finite-difference representation of the PDE hcis the conser- 
vative property, it must be established that the discretized version of the divergence 
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Figure 3.6: Eigenvalue analysis of one-dimensional model equation: fc = 1.0, A = 0.6, 
0.8, 1.0, 1.2 (periodic boundaries) 

theorem is satisfied. To do this the integral on the left is evaluated by summing 
the finite-difference representation of the PDE at all grid points. If the difference 
scheme has the conservative property, all terms will cancel except those which repre- 
sent fluxes at the boundaries. Stated another way, Eq. 3.48 holds regardless of what 
the control volume is when the ancdytical equations are modeled computationally; 
the flow field is governed by the boundary conditions. For this example, it can be 
determined whether the mass flux in equals the mass flux out. If the scheme did not 
have the conservative property, the numerical solution might permit the existence of 
mass sources and sinks [2]. 

Ni’s scheme as implemented in the present research exhibits this telescoping prop- 
erty and is indeed a conservative finite difference representation. This property is 
necessary if the flow fields to be analyzed have discontinuities such as shock waves or 
wakes. Chapter 6 investigates this idea further. 

3.6 Multiple-Grid Accelerator 

In an effort to reduce computational time associated with the solution of the linearized 
Euler equations, one would like to use a coarse grid system. The problem in doing this 
is the truncation errors described previously have been determined to be proportional 
to the square of the grid spacing. Hence, to coarse a grid leads to insufficient accuracy. 
One way to overcome this problem is to implement the multiple- grid technique [13,37] 
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Figure 3.7; Eigenvalue analysis of one-dimensional model equation: k = 0.025, 0.05, 
0.1, A = 0.8 (nonreflecting boundaries) 

where the solution on the fine grid is obtained by cycling the numerical solution 
procedure between fine and coarser grid systems. The underlying philosophy of the 
multiple-grid accelerator is to use the coarse grid to propagate the corrections of the 
fine grid at rates appreciably greater than otherwise possible on the fine grid alone. 
The result is then to obtain results with accuracy indicative of the finest grid system 
while obtaining gains in speed of convergence attributable to coarser grid structures 
without significantly adding to the computational time per iteration. 

To implement the technique, a series of coincident grids of varying spatial reso- 
lution are required. The finest grid will be designated as the zeroth level grid and 
successively coarser grids will be designated first, second, etc. level grids. These 
higher level grids are easily generated by deleting every other node in the grid on 
the previous level. Figure 3.8 provides a perspective view of three multi-grid levels. 
Recognizing that the corrections are due to wave movements and the distribution 
formulae are the systematic tools to propagate these waves, a simple multiple-grid 
scheme for solving the Euler equations is constructed by combining the basic numer- 
ical integration method with the following coarser grid solution procedure. 

Instead of the finite volume approximation as given by Eq. 3.23, the change, 
occurring in the control volumes of the 2k grid axe determined by 

= T^Huf^ (3.49) 

where T^^ is an operator which transfers to each control volume of the coarse grid the 
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Level = 2 


Level = 1 


Level = 0 


Figure 3.8; Perspective view of multiple grid acceleration levels 


correction Su’^ of the centered fine grid point, or alternatively a weighted average of 
the corrections at those fine grid points defining the coarse control volume. Making 
use of previously developed notation, the generalized distribution formulas 


(^^5U)')C' “ ^ 

Au^^ A/2^ 2 Au^'} 

(<5“nu;')c' ~ ^ 

Afl) + A^2A ^ ^y2h 

(<5«ne')c7' = \ 

Aul’i + A/“ + 


Au^') + A^2A ^ 2 


(where the prime refers to the coarse cell, C, see Fig 3.2) are then used to propagate 
the change, to the nearby coarse grid points. Next, the corrections of the 2h 

grid are found by 


= (|5«“ )^, + )c. + 


(3.51) 


After finding Au^'^ on all coarse grid points, the flow properties at the finest grid 
(level zero) are updated by 


y^new = 


( 3 . 52 ) 
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where is a linear interpolation operator which interpolates the coarse grid correc- 
tions to give the corrections at each grid point on the finest mesh. 

The above process is repeated on progressively coarser grids until the coarsest 
grid is reached. The truncation error associated with this process is second-order on 
the finest grid because the spatial approximation to the governing equation is applied 
only on the finest mesh. With multiple grid acceleration, the use of coarser grids 
reduces the computational work required to propagate unsteady disturbances out of 
the computational domain so that a steady state is rapidly reached. For a detailed 
description of the multiple grid acceleration technique the reader is referred to [13,37]. 


3.7 Summary 

In this chapter, the linearized Lax-Wendroff integration scheme used in the present 
analysis was developed. Properties of the resulting scheme, such as stability, con- 
sistency and accuracy, were also investigated. The basic Lax-Wendroff scheme was 
found to be second-order accurate and conditionally stable. In the next chapter the 
treatment of the far-field boundaries is presented. 



Chapter 4 

Far-Field Boundary Conditions 


4.1 Introduction 

Because the computational domain must be finite in extent, so-called far-field bound- 
ary conditions must be applied at the inflow and outflow boundaries. These bound- 
ary conditions can, if improperly formulated or applied, reflect unsteady disturbances 
back into the computational domain corrupting the unsteady solution. Furthermore, 
the better the far-field boundary condition, the closer the inflow and outflow bound- 
aries can be placed to the blade row of interest, thereby reducing the number of grid 
points and the corresponding computational time. 

Nonreflecting boundary conditions have been a recurring topic of investigation in 
numerical analyses, and for aeroelastic analyses in particular [19]. Previous efforts 
have focused on matching the analytical behavior of the governing equations in the 
far-field to the discretized governing equations. Verdon et. al. [49] and Whitehead 
and Grant [56] matched the known analytical far-field behavior of the linearized po- 
tential equation to finite difference and finite element representations of the linearized 
potential equation on the interior of the domain. Hall and Crawley later applied a 
similar technique to the linearized Euler equations [24]. Giles has developed approx- 
imate boundary conditions for time marching Euler applications [16]. All of these 
approaches are based on analytical descriptions of the eigenmodes of the governing 
flow equations. 

In this report, we present an alternative nonreflecting boundary condition formula- 
tion that is both more general and more accurate than previous boundary conditions 
based on analytic descriptions of the far field. The exact far-field behavior of the 
discretized equations themselves is found by solving a numerical eigenvalue problem. 
The resulting eigenmodes are then used to construct numerically exact, nonreflecting 
boundary conditions. In addition to being exact, the present formulation has the ad- 
vantage that it is generic and can be applied to other flow models (such as potential 
and Navier-Stokes equations) and can be extended to three dimensions. 

To aid in the understanding and development of the new numerically exact far-field 
boundary conditions, three separate far-field boundary conditions are presented. The 
first two are analytical in nature, i.e., the analytic behavior of the governing equations 
in the far field is determined and then coupled to the numerical integration scheme. 
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These approaches provide important insight into the far-field behavior. Finally, the 
new numerically exact far-field nonreflecting boundary conditions used in the present 
investigation is presented. 


4.2 Characteristics Of The Linearized Equations 

The following sections present the three different ways of implementing nonreflecting 
boundary conditions at the inlet and exit of the computational domain. 


4.2.1 One-Dimensional Characteristics 


The easiest nonreflecting boundary conditions to develop and apply are the so-called 
one-dimensional characteristic boundary conditions. Hence this is an appropriate 
place to start the development of the far-field boundary treatment. Recall that the 
grid motion is generated such that grid motion vanishes in the far field. This reduces 
the governing equations in the far-field to 


du' ^ ^ 

dt dx dy 


(4.1) 


If the unsteady disturbances in the far field have large cicumferential wavelengths, 
then the spatial derivative in the circumferential direction may be neglected so that 


dt dx 


(4.2) 


It will be beneficial for both the development and the understanding of the far- 
field behavior to use primitive variables rather than conservation variables in the 
subsequent analyses. To convert from the conservation variable form used thus far to 
the primitive variable form, a linear transformation is applied, i.e., 


where 


Up = [T]u' 

( P'\ 


u = 

V 


and 
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-(7-1)V 7-1 

Rewriting Eq. 4.2 in primitive variable form, and assuming that the steady flow is 
uniform in the far field gives 
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(4.4) 
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(4.6) 
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where [A] is given by 
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(4.7) 


The behavior of this system of equations (Eq. 4.6) can be determined by an eige- 
nanalysis of the matrix [A]. One can decouple the equations in Eq. 4.6 by using a 
similarity transformation. Pre-multiplying the Eq. 4.6 by the matrix of left eigenvec- 
tors, [r]~^ gives 

aw , ,aw , , 

~dt ^ ^ 

The new variables, W, in Eq. 4.8 are the so-called characteristic variables given by 


W = 


/ ^ 
W2 
W-y 

\ ^4 y 


(4.9) 


which are related to the primitive variables by 


and 


W = [Tl-'u' 


u' = [T]W 


(4.10) 


(4:ii) 


where [T] is the matrix of right eigenvectors. The diagonal matrix [A] can be written 


as 


[A] = [T]-MA][T] = 


U -a 
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(4.12) 


where a is the steady flow speed of sound. 

Each of the diagonal entries of [A] is the propagation speed of the corresponding 
characteristic wave. Consider a typical cascade operating in a regime where the axial 
Mach number is subsonic. There will be three downstream moving waves and one 
upstream moving wave. The waves have the following physical interpretation. The 
first characteristic, u;i, is given by 


W-, 


= u' — 


-pa 


(4.13) 


This corresponds to an upstream moving pressure wave. Likewise, the second chair- 
acteristic, ipj, is given by 

P' 


W2 = U' + — 

pa 


(4.14) 
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which corresponds to the downstream moving pressure wave. The third characteristic, 
1 ^ 3 , is found to be an entropy wave having the form 



(4.15) 


Finally, the fourth characteristic, W 4 , represents a vorticity wave as is given by 


W 4 = v' 


(4.16) 


Note that the vorticity and entropy waves convect with the flow, whereas the pressure 
waves propagate at the convective speed plus or minus the acoustic speed. 

Although the nonreflecting far-field boundary conditions based on one-dimensional 
characteristics can be used in some situations successfully, their effectiveness is limited 
primarily to cases with low interblade phase angles corresponding to small relative 
motions of adjacent blades. 


4.2.2 Two-Dimensional Characteristics 

Hall and Crawley [24] extended the above analysis to calculate the analytical behavior 
in the far field for the two-dimensional Euler equations. As in the one- dimensional 
approach, the grid motion in the far-field vanishes leaving the familiar unsteady per- 
turbation equation 

dn' BY' 5G' , ^ 

Agadn, it is more convenient to work with the equations in primitive variable form. 
Furthermore, it is assumed that the steady flow is uniform in the far-field so that 
Eq. 4.17 becomes 

where [A] is the matrix which appears in Eq. 4.6 and [.H] is given by 


[B\ = 


V _0 p 0 
0 V ^ 0 

0 0 y ^ 
[ 0 0 7 p y 


(4.19) 


Consider the case where the cascade is vibrating with interblade phase angle, a, 
and frequency, w. Since the solution is periodic in the circumferential direction and 
since the behavior of waves in the far-field is of interest, a more natural representation 
of the solution in the far field is given by the following Fourier series 

00 

Pm 

m=— 00 




(4.20) 
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where j = G is the blade-to-blade gap, and is a spatial wave number to 

be determined later. The coefficients of the Fourier series, are determined using 
the Fourier transform given by 



1 

— I u 
G Jo ^ 


(4.21) 


where the integral in Eq. 4.21 is evaluated numerically using the trapezoidal rule. 
Substituting Eq. 4.20 into Eq. 4.18 gives 


E 


m=— oo 


(w[/] + k„[A] + - 0 


(4.22) 


where [/] is the identity matrix and = {cr + 27rm)/G. In order for Eq. 4.22 to be 
true, each term in the series must independently vanish such that 


{o^[I]Ak^[A] + l3JB])n^^=0 (4.23) 


Since u> and are prescribed quantities, Eq. 4.23 is an eigenvalue problem for the 
eigenvalues k^ and the corresponding eigenvectors u'^. 

As in the development of the one-dimensional boundary condition, the right eigen- 
vectors can be eissembled into a matrix [T]. Then the Fourier coefficients, can 
again be written in terms of the characteristic variables, W„, . The two-dimensional 
characteristics are given by 

W„ = [T]-'u;^ (4.24) 

After some manipulation, it can be shown that 
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£a 
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(4.25) 


The matrix in Eq. 4.25 is a function of the steady flow variables, the excitation 
frequency w, and the circumferential wave numbers The resulting axial wave 
numbers are 


-U(w + + 3 + v’ - 3^) + 2^„u, V + 


nm — 




(4.26) 
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Jc*y — 


-£/(u, + l3„V-a^ 01 (U^ + V" - S’) + 20, 


,uV 


V‘-a-‘ 


(4.27) 


, _ + 

(3,4)m y 

As in the one-dimensional boundary conditions of the previous section, these charac- 
teristics correspond to upstream and downstream moving pressure waves, an entropy 
wave, and a vorticity wave. To recover the one-dimensional characteristics, a fre- 
quency, a;, of 1.0 and circumferential wave number, of 0.0 can be substituted into 
Eq. 4.25. 



4.2.3 Exact Numerical Characteristics 

Presented here is a more general method of calculating the characteristic waves for 
two-dimensional, linearized unsteady flow solvers such cis those developed for solving 
the linearized potential equation and the linearized Euler equations. Consider the 
linearized Euler equations discretized on an H-grid with I nodes in the axial direction 
and J nodes in the circumferential direction. Even though the solution procedure 
selected for this research (Chapter 3) uses a pseudo-time-marching technique and is 
iterative in nature, the converged solution can be thought of as satisfying a large 
sparse matrix equation of the form 


' B, C, 
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where u, is the vector of perturbation variables along the axial grid line, and b,- 
is the vector of inhomogeneous terms that arise from, for example, unsteady blade 
motion. As previously noted in the section detailing the deforming grid, the grid 
motion vanishes in the far-field. Hence the vectors b,- go to zero in the far field. The 
sub-matrices A,-, 5,, and C,- are large sparse matrices, each of size 4 J x 4 J. The entries 
in the these sub-matrices depend upon the details of the particular finite volume 
or finite difference scheme used in the unsteady flow solver as well as the steady 
flow solution. The details of the finite volume scheme were presented in Chapter 
3. The result of this approach is that the sub-matrices A,-, 5,-, and C, are block 
tridiagonal, with the blocks being 4x4 matrices. As previously described, upstream 
cind downstream of the blade row periodic boundary conditions are prescribed to 
reduce the computational domain to a single passage. Therefore, terms also appear 
in the upper right «ind lower left corners of the sub-matrices. 

If in the far field the grid spacing in the axial direction is uniform and the stream 
line gridlines are straight and aligned with the steady flow, then the discretized equa- 
tions are identical from axial station to axial station. The discretized equations at 



4.2. CHARACTERISTICS OF THE LINEARIZED EQUATIONS 


51 


the axial station in the far-field can be expressed a5 

[A] + [B] {uj + [C] {u.+i} = 0 (4.30) 

where now the matrices [A], [B], and [C] are independent of i. Since the equations 
are identical from station to station, a solution which satisfies the discretized far-field 
equations at the station must also satisfy the equations at the z + station. 
This, along with the insight that it is desirable to model the motion of waves in the 
far-field, suggests that solutions in the far field take on the following form 

{^i} = (4.31) 

m 

where u,- is the solution at the station, is an eigenvalue, u„i is the corresponding 
eigenvector, and is a coefficient that indicates how much of each eigenmode is 
present in the solution. Substitution of Eq. 4.31 into Eq. 4.30 gives 

[lA] + [B] + zl^ [C]] {u,„ } = 0 (4.32) 

m 

For the series in Eq. 4.32 to be zero, each term in the series must vanish so that 

[[A] + z„[B] + 410] {u„}=0 (4.33) 


This is recognized as a second order eigenvalue problem for the eigenmodes, {u„,}, 
and the corresponding eigenvalues, z„. The eigenvalue problem is put into a more 
conventional form by recasting Eq. 4.33 in state-space form. 



(4.34) 


Roughly speaking, the eigenmodes of Eq. 4.34 correspond to two-dimensional pres- 
sure, entropy, and vorticity waves which can travel up and down the duct. The 
eigenvalue is closely related to the axial wave number, of an eigenmode. In 
particular, 

zl^ = exp[ji{k^^Ax + k^rn^y)] (4.35) 


or 

, r 

Ax ’""Ax 

where Ax and Ay are the axial and circumferential shifts in the grid from one axial 
station to the next in the far-field. The matrices [A], [5], and [C] are obtained 
by careful examination of the iterative linearized Lax-Wendroff solution algorithm. 
These matrices are actually quadratic functions of the frequency of blade vibration, 
u). Rewriting Eq. 4.34 in the familiar generalized eigenvalue problem notation allows 
for a clearer presentation of the ensuing steps. Thus we write that 



[M]x = z[N]x. 


(4.37) 
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where [M] and [iV] are general complex matrices and z and x are the eigenvalue and 
right eigenvector respectively. In the present work EISPACK is used to solve this 
eigenvalue problem for the desired eigenvalues and eigenmodes. 

In figures 4.1 and 4.2 the eigenvalues computed using the present algorithm are 
compared to the analytical modes predicted using the two-dimensional analytical 
analysis presented in the previous section. Figure 4.1 shows the eigenvalues of the 
upstream far-field modes of the Lax-Wendroff scheme for a typical unsteady flow 
calculation. In this example, the blade-to-blade gap, G, is 1.0, the inflow Mach 
number, M_oo, is 0.7, the inflow angle, is 55®, the reduced frequency, k is 
1.287, and the interblade phase angle, <t, is —90°. The grid spacing between axial 
stations in the far-field in the x and y directions is 0.0601 and 0.0716, respectively. 
Figure 4.2 also presents the eigenvalues of the continuous linearized Euler equations 
as they would appear in the computational domain. Note that the eigenvalues of the 
discretized system agree well with the analytically determined eigenvalues for those 
eigenvalues near z = (1,0). These eigenvalues correspond to longer wavelength(smaller 
interblade phase angle) modes that are well modeled by the Lax-WendrofF scheme. 
Computational modes with shorter wavelengths propagate with speeds and decay 
rates different from their corresponding analytical modes. Finally in Fig. 4.1, note 
the computational eigenvalues outside the unit circle in the left half-plane. These 
modes are purely computational with no physical counterparts. 

By examining the eigenvalues of Eq. 4.37, one can determine whether the mth 
eigenmode is traveling away from or towards the rotor. For the flutter problem, 
unsteady disturbances are generated by the vibratory motion of the blade row and, 
therefore, no unsteady waves should travel toward the blade row in the far-field. At 
the upstream far-field boundary, an eigenmode with eigenvalue having a magnitude 
less than unity represents an incoming wave which decays as it progresses towards 
the blade row and hence should be excluded from the solution (this corresponds to 
the case where the wave number, k^, has a positive imaginary part). Conversely, 
if the magnitude of the eigenvalue z^^ is greater than unity, then the corresponding 
eigenmode is an outgoing mode which decays as it propagates away from the blade row 
and hence should be retained in the solution (this corresponds to the case where the 
wave number, k^^, has a negative imaginary part). If the magnitude of the eigenvalue 
z^ is unity, corresponding to the imaginary part of k^^ being zero, then additional 
analysis is required to determine the direction of wave propagation. 

In practice, eigenvalues which fall on the unit circle are one of two types: repeated 
or non-repeated roots. The repeated roots, approximately speaking, correspond to 
vorticity and entropy modes that convect with the free stream. Since the eigenvalues 
are repeated, these modes are easily identifiable as rightward moving waves (assuming 
the flow is left to right) and will be excluded upstream and retained downstream. 

The remaining modes to be analyzed have eigenvalues z^ with magnitude unity 
and are distinct. These modes correspond to upstream and downstream moving 
pressure modes that propagate unattenuated. In order to determine the direction of 
propagation of theses eigenmodes one must examine the group velocity. The group 
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Figure 4.1; Eigenvalues of upstream far-field modes of discretized equations 


velocity, is given by 

L_ 

^ dk dk/duj 

If the group velocity is positive, there is a net flux in energy to the right (incoming 
mode upstream, outgoing mode downstream). Conversely, if the group velocity is 
negative, then the net flux in energy is to the left (outgoing mode upstream, incoming 
mode downstream). 

To determine the group velocity, we first rewrite the eigenvalues as functions of 
the axial and circumferential wave numbers, 

z' = exp [j -f fcj,Ayz)] (4.39) 



Next, differentiating Eq. 4.35 with respect to w and rearranging yields 


dk^ 




1 dz 


du ■'zAxduj 

The last piece of information needed to determine the group velocity is dz/du}. To 
determine dzjdu), the following perturbation analysis is performed. Recall Eq. 4.37 
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Figure 4.2: Analytically computed eigenvalues of far-field modes 


is 

[M]x = z[iV]x 

Consider a perturbation series expansion of Eq. 4.37 about a known solution, i.e., 


([M] -I- [M']) (x -t- x') = (2 -f- z') ([iV] + [N']) (x + x') (4.41) 


where the primes indicate a small perturbation. Retaining only the first order terms 
and multiplying by the left eigenvectors, y^, gives 

{[M] - z[iV]) X' 4- y^ ([M'j - 2 [AT']) x = z'y'^[N]x (4.42) 


The first term in Eq. 4.42 is identically zero by definition. The resulting relationship 
for the perturbation of the eigenvalue, z', is 


yT{M‘ - zN')yi 
y^Nx 


Therefore, one can infer that 


dz 

du) 




y'^Nx 


(4.43) 


(4.44) 
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Finally, the group velocity can be written as 

y-l ^ +i 

3 zAx y'^Nx 

Hence, once the eigenvalues and the left and right eigenvectors of the individual 
modes are known, the group velocity and therefore the direction of propagation can 
be determined as well. 



4.3 Application of Far-Field Boundary Conditions 

Having computed the behavior of characteristic waves in the far field, we now use the 
characteristics to construct approximate and exact nonreflecting boundary conditions. 


4.3.1 One-Dimensional Nonreflecting Boundary Conditions 

First, we consider the application of the one- dimensional characteristic nonreflecting 
boundary conditions. As an example, the implementation for an upstream, far-field 
boundary node is considered. After a Lax-Wendroff iteration has been performed, 
but before the far-field boundary conditions have been applied, the estimate of the 
solution of a boundary node is given by 


u 


'n+l _ u'n ^ 
temp 


(4.46) 


Using Eq. 4.3, the solution vector in conservation form can be transformed into prim- 
itive variable form such that 


u'n-i-i 

Ptemp 



(4.47) 


In general, contains both incoming and outgoing waves. The amount of each 

of these waves is found by Eq. 4.10. With the characteristics now known, the waves 
determined to be entering the domain are zeroed since these waves could only be 
the result of a reflection at the far-field boundary. In subsonic flow, for example, 
the upstream boundary would have three characteristics entering the domain (one 
pressure, one vorticity, one entropy) which must be eliminated and one outgoing 
characteristic (pressure) that is left unchanged. This is accomplished as follows; 


= lA][ri->u 


'n+l 

Ptemp 


(4.48) 


so that 


U'n+l 


miA][T]-'u 


'n+l 

Ptemp 


(4.49) 


where [A] is a diagonal matrix with ones in the entries corresponding to outgoing 
waves and zeroes in the entries corresponding to incoming waves. Finally, to convert 
the solution back into conservation variable form 


u'n+i = 

I J Ptemp 


(4.50) 
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or 

u'n+i = (4.51) 

Downstream, the same procedure is followed except that now three characteristics 
are outgoing and must be retained and only one characteristic is incoming and need 
be deleted. 

Note that the one-dimensional characteristic nonreflecting boundary conditions 
are approximate; these boundary conditions are most effective when applied to long 
cirounferential wavelength (low interblade phase angle) disturbances. 


4.3.2 Two-Dimensional Nonreflecting Boundary Conditions 

After application of each Lax-Wendroff iteration, but before the far-field boundary 
conditions have been applied, the estimate of the solution on (for example) the inflow 
boundary is updated using Eq. 4.46. Again, the solution will contain a contri- 

bution of incoming and outgoing modes. To apply the two-dimensional characteristic 
boundary conditions, we must first determine the components of each Fourier mode 
using Eq. 4.21. Having computed the Fourier modes, two-dimensional characteristic 
nonreflecting boundary conditions are applied as follows: 


u'n+i 

Pm 


m(AllT]->u 


'n-H 

P^tcmp 


(4.52) 


where again, [A] is a diagonal matrix of zeroes and ones which eliminate incoming 
modes and [T^] the matrix of right eigenvectors. 

Once the characteristic two-dimensional boundary conditions have been applied 
to all of the Fourier modes, the modified modes are summed together using Eq. 4.20 
to obtain the solution in primitive variable form on the far-field boundary. Finally, to 
obtain the updated solution in conservation form, the following equation is applied 
at every boundary node, i.e. 

= [1]-%::^ (4.53) 

As a final note, in practice a finite number of Fourier modes are used (typically m 
= -2,-l,0,-f-l,+2). The remaining modes are not used to update the solution on the 
boundary. 

If the steady flow is uniform in the far field, then the boundary conditions de- 
scribed above are analytically exact. However, if the steady flow is nonuniform, then 
they are approximate. Also because the numerical characteristic waves differ some- 
what from the analytical characteristics due to truncation error, some small reflections 
will occur for grids with finite resolution. Finally, it is difficult to extend the two- 
dimensional characteristic nonreflecting boundary conditions to three dimensions due 
to the difficulty in obtaining analytic descriptions of three-dimensional characteristic 
waves if the mean flow contains swirl. For these reasons, a more general numerically 
exact Wcis developed and is presented in the next section. 



4.4. SUMMARY 


57 


4.3.3 Numerically Exact Nonreflecting Boundary Conditions 

In this section, we apply the numerically exact nonreflecting boundary conditions 
developed in section 4.2.3. Consider the solution along two neighboring grid lines in 
the upstream far-field region (stations i and z + 1). Together, the solution at these 
two grid lines can be thought of as the state of the solution at the grid line. In 
general, this solution will contain components of all the eigenmodes, incoming and 
outgoing. The state vector can be expressed in terms of the characteristic variables, 
W, as 

I } = [£][Aj‘W (4.54) 

where [E] is the matrix of eigenvectors found by solving Eq. 4.37 and [A] is the 
diagonal matrix of eigenvalues. Therefore the state vector at the inflow boundary is 
related to the state vector at the second axial grid line by the matrix 

where [T] is a transition matrix. The solution at the upstream boundary is then 
related to the solution on the interior by 

{-;} = [r„ (4-56) 

Finally, for gust response problems Eq. 4.56 becomes 

{«;} = K.,.} + [ 

where is a vector composed of the desired incoming eigenmodes. 

The goal is to eliminate incoming waves from the solution. Hence after each 
basic Lax-Wendroff iteration, the solution at the upstream boundary is found using 
Eq. 4.56 (for flutter calculations). The matrices Til and ^12, however, are constructed 
first by setting the entries of [A] to zero that correspond to incoming modes before 
applying Eq. 4.55. This has the effect of eliminating the incoming characteristics 
from the solution at the inflow boundary. Note that the matrices and T,, need to 
be computed just once before the start of the Lax-Wendroff iteration procedure. At 
each iteration, the far-field nonreflecting boundary conditions require only a relatively 
small matrix/ vector multiply. 

4.4 Summary 

In this chapter, three separate far-field analyses were presented culminating with 
the development of a new, numerically exax:t nonreflecting boundary condition which 
eliminates all reflections. Furthermore, although these boundary conditions were de- 
veloped for the two-dimensional linearized Euler analysis, the technique is generic and 
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can be applied to other flow models (e.g., the potential, and Navier-Stokes equations) 
and caji be to three-dimensional flows. In the next chapter, this far-field treatment, 
along with the basic numerical integration scheme presented in Chapter 3, will be 
used to investigate unsteady subsonic flows in compressors. 



Chapter 5 
Results 


This chapter presents results calculated using the present linearized Euler method. 
The cases selected show the ability of the present linearized Euler analysis to compute 
accurately and efficiently the unsteady aerodynamic response of cascade blade rows 
due to incident gusts as well as blade vibrations. The results obtained with the 
present analysis are compared with those determined using previously developed semi- 
analytical, numerical, and experimental methods. 

5.1 Flat Plate 

To validate the method, a number of unsteady flows about a cascade of flat plate air- 
foils are computed. The results are compared with those obtained using Whitehead’s 
LINSUB code [55], which is based on Smith’s compressible flat plate theory [43]. For 
all the cases considered in this section, the mean flow through the cascade is uniform 
with a Mach number, M, of 0.7. The stagger angle, 0, is 45°, and the gap-to-chord 
ratio, G, is 1.0. 

To begin, consider the case of an inlet distortion interacting with the flat plate 
cascade. In the nonrotating reference frame, the flow is axial and steady with con- 
stant total enthalpy. The axial velocity has a sinusoidally varying deficit with a 
circumferential wavelength of two blade gaps resulting in an interblade phase angle, 
<T, of —180°. In the rotating cciscade frame of reference, the airfoils see an unsteady 
gust with reduced frequency, k, (based on chord, c and upstream velocity, Vj) of 
2.221. Under these conditions the flow is superresonant, i.e., pressure waves with an 
interblade phase angle, cr, of -1-180° propagate in the far field. Any reflection of these 
pressure waves off the fcir-field boundaries would cause an unattenuated wave to prop- 
agate back into the computational domain, corrupting the solution. Therefore, this 
case provides a good test of the nonreflecting boundary conditions. Figure 5.1 shows 
the computed real and imaginary parts of the unsteady pressure difference across the 
surface of the reference airfoil. The linearized Euler results shown were computed on 
a 65 X 17 node grid and a 129 x 33 node grid. For comparison. Fig. 5.1 also shows 
the essentially exact solution calculated using LINSUB. The linearized Euler solution 
agrees very well with the exact solution, especially for the solution computed on the 
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129 X 33 node grid. This good agreement demonstrates the effectiveness of the numer- 
ically exact nonreflecting boundary conditions. Also note that the present linearized 
method accurately predicts the square root singularity of the exact solution at the 
leading edge. 

In addition to the gust response problem, one would need to calculate the aero- 
dynamic damping to complete the aerodynamic formulation of the forced response 
problem. This can be accomplished by solving the blade motion problem. Consider 
the case where the cascade of airfoils is plunging with an interblade phase angle, <r, 
of —180°, and a reduced frequency, k, of 2.221. Figure 5.2 shows that the linearized 
Euler analysis is in good agreement with the exact solution for this superresonant, 
moderate reduced frequency case. Although for the present aerodynamic damping 
analysis the computational grid deforms continuously with the specified motion of 
the airfoils, the effectiveness of the moving grid is not readily apparent since the 
mean flow is uniform and hence the error causing extrapolation terms do not appear 
in the airfoil boundary conditions. 

The last flat plate example considered is a cascade of airfoils subjected to an inlet 
distortion with an interblade phase angle, <t, of 270°, and a reduced frequency, k, of 
3.332. This corresponds to a superresonant distortion with a wavelength 1.333 times 
the blade-to-blade gap. Figure 5.3 shows the computed real and imaginary parts of the 
unsteady pressure difference across the surface of the reference airfoil. The agreement 
for this case is not as good for the previous two cases although the solution clearly 
approaches the analytical solution as the grid resolution is increased. It appears from 
these results that the scheme as currently implemented gives second-order accurate 
solutions. For example, one observes that doubling the grid resolution reduces the 
error by approximately a factor of four. Also, the errors are more pronounced at 
high reduced frequencies. This can be attributed to the fact that at higher reduced 
frequencies more resolution is required to resolve the short wavelength modes. 
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Figure 5.1; Unsteady pressure difference for fiat plate cascade subjected to an incident 
vortical gust: k = 2.221, a = —180° 



Figure 5.2: Unsteady pressure difference for flat plate cascade undergoing an unsteady 
plunging motion: k = 2.221, <r = —180° 
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Figure 5.3; Unsteady pressure difference for flat plate cascade subjected to an incident 
vortical gust: k = 3.332, a = 270° 

5.2 Tenth Standard Configuration 

Having demonstrated the accuracy of the method for a flat plate geometry, we now 
consider the unsteady flow in a more realistic compressor geometry. The purpose 
of these test cases is to demonstrate the ability of the present method to compute 
accurately the unsteady flows over loaded airfoils. In addition, the efficiency of the 
method, and the effectiveness of the far-field boundary conditions will be examined. 

The geometry considered in this section is a cascade of cambered airfoils with a 
slightly modified NACA 0006 thickness distribution. The airfoil has a circular arc 
camber distribution with a maximum height of 5 percent of the chord. The cascade 
has a stagger angle, 0, of 45° and a gap-to-chord ratio, G, of 1.0. The mean inflow 
angle, n_oo» is 55° and the inflow Mach number, is 0.7. Two computational 

grids were used for this example: a 65 x 17 and a 129 x 33 node grid. Figure 5.4 
shows the mean flow surface pressure distribution calculated using the present steady, 
nonlinear Euler solver. These results agree well with the surface pressure distribution 
computed using another steady Euler solver developed by Huff [31]. Note that the 
maximum Mach number, M, on the suction surface is about 0.92. 

Having computed the steady flow through the blade row, the unsteady flow due 
to an inlet distortion in the nonrotating frame of reference is computed. For this 
example, the interblade phase angle, <r, is —90° and the reduced frequency, k, is 1.287 
corresponding to a disturbance with a wavelength in the circumferential direction of 
four blade-to-blade gaps. In the stationary frame of reference the total enthalpy is 
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Figure 5.4: Steady coefficient of pressure distribution, Tenth Standard Configuration: 
M_«, = 0.7, G = 1.0, 0_^ = 45°, n_«, = 55° 
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constant. Figure 5.5 shows the calculated unsteady surface pressure distribution. The 
good agreement between the 65 x 17 and 129 x 33 node grid solutions indicates that 
the fine grid solution is nearly grid converged, i.e., the solution will not change with 
increased grid resolution. Note that the solution is now well behaved in the vicinity 
of the leading edge; the singularity is removed by the finite radius leading edge. Also 
shown is the solution computed using a linearized potential code (LINFLO) [26]. The 
overall agreement with LINFLO is good, although there is a slight discrepancy in the 
real part of the solution on the suction surface. 

Next, we consider a plunging motion of the compressor blades. For the first ex- 
ample, the blades vibrate with an interblade phase angle, cr, of —90®, and a reduced 
frequency, k, of 1.287. Figure 5.6 shows the unsteady pressure distribution on the 
surface of the reference airfoil found using the present linearized Euler analysis. The 
figure also includes the results of the nonlinear time marching Euler analysis of Huff. 
Figure 5.6 shows tha^th^two analyses are in excellent agreement. Note that the solu- 
tions are well behaved around the leading and trailing edges, regions that are difficult 
to resolve accurately using fixed-grid methods. These results clearly demonstrate the 
effectiveness of the deforming grid method for resolving unsteady flow features in 
regions of large mean flow gradients, even using moderate resolution grids. 

Also shown in Fig. 5.6 are linearized Euler solutions calculated on fixed (i.e., 
nondeforming) computational grids. For these solutions, extrapolation terms which 
depend on mean flow gradients must be added to the airfoil boundary conditions and 
to the expression for the unsteady surface pressure to account for the fact that the 
airfoil vibrates through the stationary grid. Large errors are seen in the solutions 
computed on the fixed grid, especially around the leading and trailing edges. These 
errors are inevitable when using a fixed grid due to the difficulty in evaluating the 
gradient of the mean flow field near the airfoil surface. 

Finally, we consider a moderately high reduced frequency blade motion case to 
demonstrate the current limitations of the method. For this case, the airfoils plunge 
with an interblade phase angle, a, of —180°, and a reduced frequency, k, of 2.573. 
Figure 5.7 shows the computed unsteady surface pressure distribution. Also shown for 
comparison are the results of a linearized potential analysis [21] and a time marching 
Euler calculation. The potential solution, which wais computed on an extremely fine 
gird (200 x 50 nodes), is grid converged. The fine grid linearized Euler calculation 
is in very good agreement with the potential calculation. The coarse grid linearized 
Euler calculation, on the other hand, differs significantly from the potential solution 
on the suction surface. This is to be expected since on the suction surface, where 
the Mach number is large, upstream travelling pressure disturbances have very short 
wavelengths that axe difficult to resolve. Of course, adequate grid resolution is essen- 
tial for both nonlinear and linearized analyses. 

Figure 5.8 shows the computed unsteady pressure contours for this blade motion 
case. Of particular interest is the behavior of the solution in the feir-field region. For 
this example, the flow is superresonant upstream and subresonant downstream. In 
both the upstream and downstream far-field regions, the pressure contours aire seen 
to pass smoothly out of the domain without reflection demonstrating the effectiveness 
of the far-field nonreflecting boundary conditions. 
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Figure 5.5: (top) Real unsteady pressure distribution for cascade of Tenth Stan- 
dard Configuration airfoils subjected to an inlet vortical gust; (bottom) Imaginary 
unsteady pressure distribution for cascade of Tenth Standard Configuration airfoils 
subjected to an inlet vortical gust: k = 1.287, cr = —90° 
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Figure 5.6: (top) Real unsteady pressure distribution for cascade of Tenth Standard 
Configuration airfoils undergoing an unsteady plunging motion; (bottom) Imaginary 
unsteady pressure distribution for cascade of Tenth Standard Configuration airfoils 
undergoing an unsteady plunging motion: k = 1.287, cr = —90° 
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Figure 5.7: (top) Real unsteady pressure distribution for cascade of Tenth Standaxd 
Configuration airfoils undergoing aji unsteady plunging motion; (bottom) Imaginary 
unsteady pressure distribution for cascade of Tenth Standard Configuration airfoils 
undergoing an unsteady plunging motion: k = 2.573, a = —180® 
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Figtire 5.8: Unsteady pressure contours for cascade of Tenth Standard Configuration 
eurfoils undergoing an unsteady plunging motion: k = 2.573, a = —180" 
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To further substantiate the effectiveness of the far-field boundary conditions, the 
same case was computed again. This time, however, the far-held boundaries were 
placed farther away from the blade row of interest. The unsteady surface pressures 
computed in the lengthened computational domain are compared to those previously 
computed and are shown in Fig. 5.9. The results indicate that the close proximity of 
the far-field boundary does not affect the unsteady pressure distribution. Addition- 
ally, by moving the boundaries closer to the blade row of interest, fewer computational 
nodes and hence computational resources are required to solve the problem. 

Figure 5.10 shows the convergence histories for this flutter calculation (on the 
deforming grid) using the basic Lax-Wendroff solver alone and with Ni’s multiple- 
grid acceleration technique. The convergence criteria is that the L2 norm of the 
pu corrections be less than 5.0 x 10~®. Without the multiple-grid accelerator, the 
linearized unsteady Euler analysis required 6433 iterations to converge. With three 
levels of the multiple-grid accelerator, the linearized unsteady Euler analysis required 
just 773 iterations corresponding to about 14 minutes on a Stardent 3000 workstation. 
The nonlinear time marching solution, on the other hand, required 7384 iterations 
and approximately 88 minutes on a CRAY-YMP. When taking into consideration the 
difference in speeds of the two computers, one concludes that the linearized Euler 
analysis is nearly two orders-of-magnitude faster than the nonlinear analysis while 
still modeling the essential features of the unsteady flow. Additionally, for cases 
considered thus far, the multiple-grid accelerator works as well for the linearized 
unsteady Euler analysis as for the nonlinear, steady Euler analysis. The reduction 
seen in iteration count by additional multiple grid levels is approximately the same 
for both the steady and unsteady analyses, and both require approximately the same 
number of iterations to reach convergence. 

Finally, having computed the unsteady pressure distribution on the airfoil surface, 
one can integrate to obtain the unsteady pitching moment and hence draw a conclu- 
sion about the flutter stability of the cascade (a positive imaginary part corresponds 
to negative cierodynamic damping). Shown in Fig. 5.11 is the imaginary part of the 
pitching moment for a cascade of Tenth Standard Configuration airfoils computed 
for a range of interblade phase angles, cr, from —180° to -|-180°, and at a reduced 
frequency, k, = 0.5. These results were obtained using a 65 x 17 node computational 
grid. Also shown is the pitching moment computed using Hall’s linearized potential 
method. The peaks in the solution are acoustic resonances. Note that there is gen- 
eradly good agreement between the linearized Euler and linearized potential solution. 
The good agreement (particularly in the superresonant regions) further substantiates 
the effectiveness of the far-field boundary conditions. 

5.3 First Standard Configuration 

In order to further validate the present method, we consider the Cctse of a low speed 
compressor cascade known as the First Standard Configuration. This cascade was 
studied experimentally by Carta [10] and provides a good test of the present method’s 
ability to predict the unsteady aierodynamic loading due to blade motion. The cascade 
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Figure 5.9: (top) Real unsteady pressure distribution for cascade of Tenth Standard 
Configuration airfoils undergoing an unsteady plunging motion; (bottom) Imaginary 
unsteady pressure distribution for cascade of Tenth Standard Configuration airfoils 
undergoing an unsteady plunging motion: k = 2.573, a = —180° 
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Figure 5.10: Convergence histories for cascade of Tenth Standard Configuration air- 
foils undergoing an unsteady plunging motion: k = 2.573, <r = —180° 
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Figure 5.11: Imaginary part of unsteady pitching moment for cascade of Tenth Stan- 
dard Configuration airfoils undergoing an unsteady pitching motion: k = 0.5 

is composed of airfoils which have 10® of camber and are 6 percent thick. The stagger 
angle, 0, is 55® and the gap-to-chord ratio, G, is 0.75. The reported inflow angle, 
fl_oo 5 is 66® and the inflow Mach number, M_^ , is 0.17. Figure 5.12 shows the 
steady surface pressure distribution calculated with the present nonlinear, steady 
Euler solver in addition to the experimentally measured values. For the unsteady 
case considered here, the airfoils vibrate in pitch about a point near their midchords 
with a reduced frequency, k, of 0.244 and an interblade phase angle, <7, of —90°. The 
computed unsteady surface pressure is shown in Fig. 5.13 along with the experimental 
results of Carta. The computational results were obtained using a 129 x 33 node grid. 
Overall, the agreement between the present analysis and the experiment, while not 
exact, is qualitatively very good. 

5.4 Summary 

The results obtained using the present linearized Euler analysis validated the present 
method agciinst previous analytical, numerical, and experimental results. First, the 
flat plate cases studied were compared with an analytical solution computed using 
Whitehead’s LINSUB. The excellent agreement proved that at least for flat plate cas- 
cades, the present linearized analysis can accurately and effectively calculate the un- 
steady aerodynamic forces for both the forced response and flutter problems. Next, a 
cascade of loaded airfoils (the Tenth Standard Configuration) was analyzed for several 




Steady Coefficient of Pressure, Cp 
0.80 0.40 0.00 -0.40 -0,80 


5.4. SUMMARY 


73 




Distance Along Chord, ^/c 


Figure 5.12: Steady coefficient of pressure distribution, First Standard Configuration: 
M_^ = 0-17, G = 0.75, = 55°, = 66° 
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Figure 5.13; (top) Real unsteady pressure distribution for cascade of First Standard 
Configuration airfoils undergoing an unsteady pitching motion; (bottom) Imaginary 
unsteady pressure distribution for cascade of First Standard Configuration airfoils 
undergoing em unsteady pitching motion: (case 8) k = 0.244, <7 = —90° 
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unsteady operating conditions. The results agreed well with Hall’s linearized poten- 
tial as well as Huff’s time-accurate time-marching unsteady Euler analysis. Finally, 
results from the present method were compared with Carta’s experimental data (First 
Standard Configuration). Based on these comparisons, the present analysis appears 
to be accurate for a variety of geometries and flow conditions. In addition to accu- 
racy, the results accompanying the Tenth Standard Configuration illustrated several 
important features of the present method. These include: (1) the far-field boundary 
conditions effectively eliminate spurious reflections, (2) the continuously deforming 
unsteady grid has effectively eliminated error producing extrapolation terms inherent 
in fixed grid calculations, and 3) the present linearized analysis has taken advantage 
of traditional computational acceleration techniques such as multigrid acceleration 
and local time-stepping. The conclusions drawn from these results are as follows. 
First, linearized anailyses are indeed a viable method for the accurate prediction of 
unsteady loading in turbomachinery applications. Secondly, the computational effi- 
ciency makes linearized methods preferable to time-accurate time-marching methods 
as long as the sm^lll motion assumption of the linearized analyses is not violated. 



Chapter 6 

Transonic Theory 


6.1 Introduction 

In recent years, a number of linearized flow analyses have been developed to compute 
unsteady flows in cascades, especially the unsteady flows that produce the aeroelastic 
phenomena of flutter and forced response. The unsteady aerodynamic loads acting 
on trainsonic airfoils in cascades are composed of two parts: the unsteady pressure 
distribution away from the shock, and a “shock impulse” load that acts where the 
shock impinges on the airfoil surface. This shock impulse arises from the unsteady 
motion of the shock. Accurate prediction of the shock impulse is important since 
the unsteady aerodynamic load due to the shock impulse is of the same order as the 
unsteady aerodynamic loads due to the unsteady pressure away from the shock. In 
viscous flows, the shock is smeared near the airfoil surface due to shock/boundary 
layer interaction and hence, strictly speaking, no shock impulse exists at the surface. 
Away from the airfoil, however, the shock wave is very thin, typically on the order 
of a few mean free paths thick, and the concept of a shock impulse is important in 
connecting the regions of smooth flow on either side of the shock. 

Verdon et al [48,50] and Whitehead [51,55] have developed linearized potential 
analyses of two-dimensional subsonic and transonic flows in cascades. Both Verdon 
and Whitehead have used shock capturing to model unsteady shock loads. Verdon 
has also used shock fitting in his linearized potential analysis to explicitly model the 
shock motion. Because of the assumption of isentropic and irrotational flow, however, 
these potential analyses cannot be used to model unsteady flows with strong shocks, 
flows with shocks that span the blade passage, or general three-dimensional flows. For 
this reason, investigators have begun to develop linearized Euler analyses of unsteady 
cascade flows [22,24,29,32]. Hall and Crawley [24] have shown that shock fitting can be 
implemented within the framework of a linearized Euler analysis to model accurately 
the unsteady motion of shocks. However, due to the inherent complexity of shock 
fitting algorithms, one would prefer to use the simpler shock capturing technique to 
model the shock impulse. 

While shock capturing is favored for its simplicity, it has only recently been shown 
that the shock impulse load can be modelled properly using shock capturing within a 
linearized framework. There are two approaches that have been suggested for obtain- 
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ing discretizations of the linearized Euler equations. The first approach, referred to in 
this chapter as Method I, is to first discretize the nonlinear unsteady Euler equations 
and then linearized the resulting finite difference equations. The second approach, 
Method II, is to first linearize the nonlinear unsteady Euler equations, then discretize 
the resulting linearized equations using traditional finite difference or finite volume 
techniques. Lindquist and Giles [35,36] have argued that the unsteady shock loads 
will be correctly predicted provided the linearized code is a Method I type lineariza- 
tion of a time-accurate, conservative, nonlinear flow solver. Their results thus far, 
however, have been limited to quasi one- dimensional channel flows. Furthermore, 
they do not discuss the conditions under which Method II linearizations will properly 
model the shock impulse. 

The objectives of this chapter are twofold. First, we demonstrate mathematically 
and by numerical experiment that the requirement put forth by Lindquist and Giles 
that the linearization be a Method I linearization of an unsteady nonlinear scheme 
is too stringent. We show that Method II linearizations will also work so long as 
the finite difference representation of the linearized Euler equations is conservative. 
Second, having demonstrated that conservative Method II linearizations may be used 
to properly model the unsteady shock impulse, we present a linearized Euler analysis 
(Method II type) of unsteady two-dimensional flow in cascades. Ni’s Lax-Wendroff 
scheme [37] is used to obtain a finite volume representation of the unsteady linearized 
Euler equations. Computational results are presented for both and two- and three- 
dimensional unsteady transonic flows in cascades. Some of these calculations are 
compared to those computed using a nonlinear time-marching shock capturing Euler 
analysis. It is shown that the present unsteady linearized analysis agrees quite well 
with the nonlinear analysis, and further that the present linearized analysis is nearly 
two orders-of-magnitude more efficient than the nonlinear analysis. The computed 
results also demonstrate that the unsteady shock loads can provide a destabilizing 
influence on the flutter stability of cascades. 

6.2 Theory 

6.2.1 Flow Field Description 

In this chapter, we again assume that the unsteady flow is inviscid and adiabatic, and 
that the unsteady flow in a cascade may be modelled by the Euler equations. For a 
two-dimensional Cartesian coordinate system, the Euler equations are given by 

^ ^ ^ 

dt ^ dx dy 


( 6 . 1 ) 
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A, A A 

where U is the vector of conservation variables, F, and G are the so-called flux 
vectors. These vector quantities are given by 
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pu 


pv 

pu 

, F = 

piP -f- p 
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puv 

pv 


puv 


ptj2 -1- p 

e 
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t-C 


where p is the density, p is the pressure, u and v are the x and y components of 
velocity, is the internal energy, and h is the enthalpy. The pressure, p, and the 
enthalpy, h, are given by 


p = (7 - 1) 


e - -p + i)2) 


and 

h = 

Next, we would like to determine the smedl disturbance behavior of Eq. 6.1 due 
to, for example, the fluttering motion of the blades of the cascade. To improve the 
accuracy of these calculations, a number of investigators have proposed the use of 
a harmonically deforming computational grid [22,23,29]. The motion of the grid is 
defined by 

4 + /(6 (6-2a) 

= 7 + (6-2b) 

^ (6.2c) 


e + P 7 P . 1 

= TT + o (« + ^ ) 

P 7 - 1 P 2 


where w is the frequency of vibration of the blades, and where / and g are the 
perturbation amplitudes of the grid motion about the mean positions, ^ and p. Having 
defined the grid motion, the unsteady flow field is represented by the perturbation 
series 

U(^, 77 , r) = U(^, 7j) -1- u(^, (6.3) 

Substitution of Eq. 6.3 into Eq. (6.1) and collection of the terms that are first-order 
in the perturbation u results in the linearized Euler equations, 

d fdF \ d /dG \ ^ 

where b is a fairly complex expression which depends on the mean flow and the 
prescribed grid motion. 


6.2.2 Numerical Modelling of the Shock Impulse 

The first question we address in this chapter is; What is the proper way to discretize 
and linearize the Euler equations in such a way that the linearized finite difference or 
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finite volume equations properly predict the shock impulse loads that result from the 
unsteady shock motion. There are two obvious approaches one can take to obtain a 
discretization of the linearized Euler equations. One approach (Method I) is to first 
discretize the nonlinear Euler equations and then linearize the resulting nonlinear 
finite difference equations. The other approach (Method II) is to first linearize the 
nonlinear Euler equations then discretize the resulting linearized equations. We claim 
here that both approaches will produce the correct result provided that the resulting 
difference equations are conservative (we give a more precise definition of what con- 
servative means in the linearized case shortly). A mathematical justification of this 
conjecture is given below. 

Due to the complexity of the two-dimensional Euler equations, we consider the 
simpler one-dimensional model equation given by 


dU dF „ 

"a + 


(6.5) 


where F = F{U), P = P{tl), and B = B{x). This model equation is very similar in 
form to the quasi-one- dimensional Euler equations which describe flow in a channel 
with a spatially varying cross sectional area. Since F and P are in general nonlinear 
functions of the conservation variable U, this model equation is nonlinear. 

As before, we model the conservation variable U as the sum of a mean part U 
plus a small harmonic perturbation ue-'"*. The mean solution is solution is governed 
by 

dF 

— + BP = 0 (6.6) 

where F = F{U) and P = P{U). The linearized unsteady model equation is given by 


d fdF 


JUJU + —\ 


dx \dU 


u 




(6.7) 


where u is the perturbation solution, and and are steady flow Jacobians. 

Returning for the moment to the unsteady nonlinear model equation, Eq. 6.5, it 
is well known that because the model equation is nonlinear, it will in general admit 
genuine solutions, that is, solutions with flow discontinuities. In smooth regions of 
the flow, the genuine solutions satisfy the differential equation, Eq. 6.5. The weak 
solution is that genuine solution which also satisfies the integral relation 


J j {ot^ ~ 9xP^ J g{x,0)U{x,0)dx = 0 (6.8) 


for every test function g(x,t) which vanishes for large x or f and which has continuous 
first derivatives [33,34]. One can then show that the unsteady Rankine-Hugoniot 
shock jump conditions at flow discontinuities are given by 




(6.9). 
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where the symbol [[• • •]] denotes the jump in the enclosed quantity across the shock, 
and Xg is the velocity of the shock. 

If one then considers an unsteady flow with small harmonic unsteadiness, one may 
linearize Eq. 6.9 to obtain the linearized shock jump conditions [24]. 


jujx, [[17]] - 


\\dF 

U 

— 

■5F1' 

[[dU 


_ 5x J 


= 0 


( 6 . 10 ) 


where u is the small disturbance part of the unsteady flow and is the small complex 
amplitude of the shock motion. Noting that the steady flow solution is given by 
dFjdx = —BP, Eq. 6.10 may be rewritten as 

'dF 


jux, [[{/]] - 


dU 


u 


+ x,B [[Pll = 0 


(6.11) 


A graphical interpretation of Eq. 6.11 is shown in Fig. 6.2.2. Shown are the mean 
and unsteady flow shock trajectories as well a.s the resulting unsteady flow, O’, the 
mean flow U, and the perturbation flow, u. The latter is just the difference between 
the unsteady and mean flows, u = 0 — U. Note that near the shock, an impulse in u 
appears due to the motion of the shock. In the limit as the unsteadiness in the flow 
tends toward zero, the integrated value of this impulse is given by 

rXs+e 

/„ = / udx = -r, [[t/]] = -X, (U^ - U,) (6.12) 

^ Xc— € 


fXs~c 

Finally, Eq. 6.11 may be written as 

jujiu + 


IdF 

du 


u 


+ B/, = 0 


(6.13) 


We presently demonstrate that the weak solution of the linearized unsteady model 
equation, Eq. 6.7, produces an equivalent shock jump condition. Multiplying Eq. 6.7 
by a test function ^(x) and integrating the result over the solution domain x 6 [0,T], 
we obtain 

u ) + B—u\ 


I 


9{x) 


juju + ^ h 


dx = 0 


dx \dU~J ' ~dU 
Integration by parts applied to the middle term in Eq. 6.14 gives 


I 


^0 dgdF 'dP 
J^gu-—^^u + gB—u 




= 0 


(6.14) 


(6.15) 


dxdU" ' ^ dU 

Next we let the test function 5 r(x) be given by 

«<*> = { J o 

where X, is the mean shock position and e is a small positive number. Differentiating 
Eq. 6.16 with respect to x gives 


if Xj - € < X < X, + e 

otherwise 


(6.16) 


^ = 6[i - {X. - e)l - S[x - {X, + £)] 


(6.17) 
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where is the Dirac delta function. Substitution of Eqs. 6.16 and 6.17 into 

Eq. 6.15 gives the desired shock jump conditions of the linearized unsteady model 
equation, 


rX,+t 


ju 


udx + 


IXs-l 


dF 

dU 
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+ B 


rXs+c Qp 


udx = 0 


(6.18) 


The integrals in Eq. 6.18 are the areas under the impulse in u and the impulse in 
p, the perturbations in the conservation variable and pressure, respectively, and are 
denoted here by and 7p (See Fig. 6.2.2). 

Finally then, we may write the Rankine-Hugoniot jump conditions for the lin- 
earized unsteady Euler equations as 


+ 


dF 

W 


u 


+ BIp — 0 


(6.19) 


This expression is identical to Eq. (13) thus demonstrating that the weak solution to 
the linearized unsteady model equation is the same as the linearized weak solution of 
the nonlinear model equation. 

We conclude, therefore, that for a finite difference scheme to properly model the 
linearized unsteady model equation, the finite difference scheme must be stable and 
consistent and satisfy the condition given by Eq. 6.15 in the limit as Ax and At tend 
toward zero. In other words, the order of linearization is immaterial; what matters 
is whether the resulting discretization is conservative. Note that this condition is 
less stringent that the condition suggested by Lindquist and Giles [35,36] that the 
discretization be both conservative and a Method I linearization. Also note the im- 
portance of the area of the impulse. When capturing the shock impulse, the width 
and height of the impulse will depend on the amount of smoothing (or artificial vis- 
cosity) in the numerical scheme. The area under the impulse, however, should be 
independent of the smoothing. 


6.2.3 Method I and Method II Linearizations 

To illustrate the difference between Method I and Method II linearizations, we con- 
sider again the model equation given by Eq. (5) with the source term set to zero 

{B = 0). 

Method I Linearization 

Consider the discretization of the nonlinear unsteady model equation, Eq. (5), using 
the Lax-Wendroff scheme. The one-dimensional computational grid is assumed to 
have constant cell size Ax and constant time step At. The solution at time level 
n -h 1 is found by Taylor expanding the solution about time level n to obtain 

-f 0{At^) (6.20) 


^n+l ^ IJn ^ 


dt 
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At2 d^U\ 
2 dt^ 
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Figure 6.1: Top to bottom: (a) Trajectory of shock in a channel or on airfoil sur- 
face; (b) Mean and unsteady flow distribution; (c) Perturbation flow showing shock 
impulse; (d) same as (c) with impulse modelled by shock capturing. Note the area 
under the impulse is the same in (c) and (d). 
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where i denotes the ith grid node in the x-direction. The time derivatives in Eq. 6.20 
are obtained by manipulation of the original model equation, Eq. (5). Rearranging 
Eq. (5) gives 


dU _ dF{U) 

dt dx 

Differentiating Eq. 6.21 with respect to time gives 


( 6 . 21 ) 


di^ dx di j dx ydU dx J 


( 6 . 22 ) 


Next, substitution of Eqs. 6.21 and 6.22 into Eq. 6.20 yields 


« , . dF 

jyn+l = [/n_ 
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(6.23) 


Finally, using centered finite difference expressions to approximate the spatial deriva- 
tives, the familiar Lax-Wendroff finite difference equation is obtained, i.e.. 


^n+l ^ Ifn, 
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^ t-Hl «-l 


AP 
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(6.24) 

The Lax-Wendroff scheme above is second-order accurate [0{Ax'^, Ai^)] and is con- 
servative (since by assumption Ax and At are constant). A Fourier stability analysis 
indicates the scheme is stable for CFL numbers less than unity. Lax [33] and Lax 
and Wendroff [34] have shown that if the conservative form of the Euler equations is 
used, and the discretization of the Euler equations satisfies numerical conservation, 
and further that the scheme is consistent and stable, then the shock wave speed and 
strength will be correctly predicted. The Lax-Wendroff scheme above satisfies these 
conditions and therefore will correctly predict the unsteady shock motion. 

Next, we would like to determine the behavior of Eq. 6.24 when the unsteadiness 
in the flow is small compared to the mean flow field. Since the resulting equations will 
be linear, we may without loss of generality assume that the flow field is composed 
of a nonlinear steady mean flow and a small perturbation harmonic unsteady flow so 


that 

U{x,t) = U{x) + u{x)e^‘^ (6.25) 


where u is much smaller than U. When viewed on our computational grid, Eq. 6.25 
becomes 

Ur^ = Ui -t- (6.26) 

Substitution of Eq.6.26 into Eq. 6.24 and collecting terms of first order in the per- 
turbation quantity u gives the desired discrete small disturbance behavior of the 
nonlinear finite difference equations. 
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(6.27) 


Equation 6.27 describes the small disturbance behavior of the nonlinear Lax-WendrofF 
equation. One interesting feature of Eq. 6.27 is the appearance of the terms involving 
d'^FIdU'^. These terms appear because of nonlinearities of the Lax-Wendroff scheme 
itself rather than nonlinearities of the Euler equations. 

For a one-dimensional problem, Eq. 6.27, along with appropriate inflow and out- 
flow boundary conditions, could be assembled into a tridiagonal matrix equation 
which could then be solved quite efficiently using Gaussian elimination (the Thomas 
algorithm). For two- and three-dimensional problems, however, this approach would 
be computationally expensive and require large amounts of computer storage. For 
these reasons, an iterative solution technique is preferred. The following explicit 
relaxation procedure is proposed: 


„n+l = q. (0 28) 

where is the left-hand side of Eq. 6.27. As Eq. 6.28 is marched in time, a steady 
state value of uj* will be obtained and the solution to Eq. 6.27 will be recovered. 
This procedure is similar to the pseudo-time time-marching technique proposed by 
Ni and Sisto [39] for solving the linearized Euler equations. Equation 6.27 can be 
shown to be consistent with the linearized model equation, Eq. (7), with truncation 
errors which are C?( Ai^, At^). A Fourier analysis of Eq. 6.28 reveals that the scheme is 
unconditionally unstable if is non zero. A spectral radius stability analysis, however, 
that taJces into account the stabilizing effect of the far-field boundary conditions shows 
that the scheme is stable for CFL numbers less than unity [12]. 

Method II Linearization 

An alternative approach to Method I is to first linecirize the nonlinear unsteady flow 
equations, an d the n discretize^ the resulting linear equations. To illustrate this ap- 
proach, we return again to the one-dimensional model equation given by Eq. (5) 
and, introducing the pseudo-time assumption of Ni and Sisto [39], assume that the 
unsteady flow U{x,t) is composed of a nonlinear mean flow, U{x), plus a small un- 
steady harmonic perturbation flow, u{x, so that 

U{x, t) = U(x) + u{x, t)eJ-‘ (6.29) 

Substitution of Eq. 6.29 into Eq. (5) and collection of first-order terms results in the 
pseudo-time linearized model equation 
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Note that Eq. 6.30 is now hyperbolic in time so that it can be marched in time. 
Furthermore, as time advances, u will reach a steady-state value so that the solution 
to Eq. (7) will be recovered. 

The next step is to discretize the linearized equation using the Lax-WendrofF 
scheme. Manipulation of Eq. 6.30 gives 


du d {dF \ 

_ = _ _ l—uj 


(6.31) 
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Finally, making use of centered spatial derivatives and substitution into the Taylor 
expansion, Eq. (20) gives the desired Lax-WendrofF formula, 
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(6.33) 


As in the Method I discretization [Eq. (27)], the Method II discretization [Eq. 6.33] is 
consistent with the linearized model equation [Eq. (7)] with truncation errors which 
are 0{Ax^, At"^). 

Note that the Method II discretization [Eq. 6.33] differs significantly from that 
obtained using Method I [Eq. (27)]. In particular, the unsteady terms involving u 
are somewhat different, and the quasi-steady terms involving d'^FjdU'^ in Eq. (27) 
do not appear in Eq. 6.33. Clearly, the order in which the linearization is performed 
is important in determining the precise form of the difference equations. 
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6.2.4 Test for Linearized Conservation 

Consider the Method I discretization of the model equation, Eq. (27). The one- 
dimensional computational grid has M nodes. To test for conservation, Eq. (27) 
is multiplied by p,Ai/At (where = gi^i)) and summed over the computational 
domain. After some manipulation including summation by parts, one can show that 
this sum is given by 
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(6.34) 

In the limit as At, Ax 0, Eq. 6.34 approaches Eq. (15). Therefore, this Method 
I linearization of the Lax-Wendroff scheme is conservative (at least for the case con- 
sidered here of constant At, Ax). A similar analysis of Method II reveals that is also 
conservative. Hence, both methods should correctly predict the shock impulse. 


6.3 Results 

6.3.1 Transonic Channel Flow 

To test the present linearized Euler analyses, we first consider the transonic flow 
through a diverging channel. This case is presented to demonstrate the ability of 
the linearized Euler method to model shock motion accurately using shock capturing. 
We will demonstrate that both Method I and Method II linearizations will produce 
satisfactory results so long as they are conservative. 

The channel considered here has a height, A, given by, 

A{x) = Ai„,^t 1 1.10313 -1- 0.10313 tanh 

(The units may be taken to be any consistent set of units.) So that we may compare 
the results obtained by the present method to those obtained by a one-dimensional 
shock-fitting theory, A^^ut is taken to be small compared with the channel length 
(■^inUt = 0.01). The inflow total pressure, Pj, total density, px, and flow velocity, 
U, are 1.0, 1.364, and 1.0 respectively. The back pressure, Pexiti is 0.7422. Shown in 
Fig. 6.3.1 is the Mach number and pressure distribution as computed using the present 
nonlinear steady Euler solver on a 129 x 5 node computational grid. The grid was 
generated so that the computational cells all have the same area, AA. The time step. 
At, used in these calculations weis constant throughout the computational domain 
unless otherwise noted. Constant Ai and AA were chosen because Ni’s scheme is only 
conservative if the ratio AtfAA is constant throughout the computational domain. 

Also shown for comparison in Fig. 6.3.1 is the solution determined using a steady 
quasi-one-dimensional, shock-fitting Euler solver using 1001 grid nodes in the i- 
direction. The shock-fitting Euler solution is grid converged and may be taken to 


10(x-- 


,0 < z < 1 


(6.35) 
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be the exact solution. Note the excellent agreement between the two different ap- 
proaches. The only noticeable differences occur at the shock, where the present 
nonlinear Euler analysis distributes the shock over about five grid nodes. 

Next, we consider a quasi-steady perturbation in the back pressure. The pertur- 
bation solution was calculated using four different approaches. First, the solution 
wais calculated using a quasi-one-dimensional, shock-fitting, linearized Euler analysis. 
This solution was computed on an extremely fine grid (1001 nodes in the x-direction) 
and is essentially the exact solution. Next, the present nonlinear steady Euler solver 
was used to compute two nonlinear solutions at slightly different back pressures. 
These two solutions were then subtracted one from the other and the result was 
normalized by the difference in back pressures to obtain the perturbation solution. 
Finally, the solution was determined using the present linearized Euler analysis (both 
Methods I and II). It should be noted that for this comparison the usual nonreflecting 
boundary conditions were replaced with reflecting boundary conditions. Upstream the 
perturbation in total pressure and density as well as the inflow angle was set to zero. 
Downstream the perturbation in static pressure was prescribed. These boundary con- 
ditions for this model problem were chosen for their simplicity and are not meant to 
model any real physical system. The results of these various approaches are shown in 
Fig. 6.3.1. As expected, all of the solutions are in excellent agreement in regions away 
from the shock. At the shock, however, the methods using shock-capturing produce 
an impulse of pressure. This impulse arises from the shock displacement. The area 
under the impulse is equal to the product of the shock displacement and the mean 
pressure jump across the shock. The shock impulse then represents the load exerted 
on the wall due to the motion of the shock. Also shown in Fig. 6.3.1 is an enlarged 
view of the shock region. Note that the computed results from the Method I and 
Method II linearizations are virtually identical to the perturbation of the nonlinear 
Euler analysis. 

To further validate the linearized shock capturing technique for unsteady flows, we 
computed the unsteady pressure distribution due to an unsteady perturbation in back 
pressure with an excitation frequency, w, of 1.0. The results are shown in Fig. 6.3.1. 
The Method I and Method II results are indistinguishable from one another and are 
therefore plotted with a single symbol. Also shown are the results of a quasi one- 
dimensional, unsteady, shock-fitting, linearized Euler solver. Away from the shock, 
the results agree quite well with the Method I and II results. At the shock, the present 
Method I and II solutions show an impulse. This impulse represents the unsteady 
load acting on the channel wall due to the motion of the shock. 

To determine whether the present linearized Euler solver correctly predicts the 
unsteady loads induced by the shock motion, the pressure was integrated over the 
lower channel wall to determine the net wall force. The results from this anedysis are 
tabulated in Table 6.1 for several different frequencies. Also tabulated in Table 6.1 
is the wall force computed using the linearized unsteady shock-fitting code. The 
agreement between the conservative form of the Method I and Method II analyses are 
seen to be in almost perfect agreement with the shock fitting scheme for all frequencies 
suggesting that the shock impulse found using shock capturing is properly modelled. 
Even in the higher frequency cases the agreement is quite good, although there is 
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Figure 6.3: Top: Perturbation pressure in a diverging channel due to a steady per- 
turbation in back pressure. Bottom: Enhancement of the shock impulse region. 
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Figure 6.4: Top; Unsteady pressure in a diverging channel due to an unsteady per- 
turbation in back pressure, u = 1.0. Bottom: Enhancement of the shock impulse 
region. 
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Table 6.1: Predicted pressure loads in a transonic diverging channel due to an un- 
steady perturbation in back pressure using a uniform axea computational grid. 


Frequency, u 

Scheme 

Wall Force 


0.0 

ID Shock Fitting 
Nonlinear Euler® 
Method I 
Method II 

1.0305 Z 0.0® 
1.0346 Z 0.0° 
1.0273 Z 0.0® 
1.0273 Z 0.0° 


1.0 

ID Shock Fitting 
Method I 
Method II 

0.6390 Z -78.7® 
0.6353 Z -78.8® 
0.6354 Z -78.8^ 



Method I 

0.5397 Z -46.6° 

Noriconservative ** 


Method II 

0.6229 Z -84.9° 

Nonconservative ® 

2.0 

ID Shock Fitting 
Method I 
Method II 

0.1974 Z - 114.1° 
0.1983 Z -113.6° 
0.1984 Z - 113.6° 



“Results from the steady analysis were found for two slightly different back pressures. 
The two solutions were then differenced and normalized by 

^Time accurate time marching steady and unsteady solution on a nonuniform area 
computational grid. 

“Local time stepping used in steady and unsteady analyses. 


a slight error (about 0.5®) in the phase of the wall force. It is believed that these 
dilferences arise from the dispersion errors in the solution away from the shock rather 
than from a limitation in shock capturing at high reduced frequencies. 

Finally, for the u; = 1.0 case, we deliberately made the Method I and II calculations 
nonconservative to demonstrate that the shock impulse cannot be properly modelled 
using a nonconservative algorithm. In the Method I calculation, the time step At 
was held constant throughout the domain, but a grid with variable cell areas near 
the shock was used. In the Method II calculations, the cell areas Aj 4 were constant 
throughout the computational domain, but the time step used in each computational 
cell was based on a local CFL number (local time stepping). In both cases, the ratio 
AtfAA varies over the computational domain making the schemes nonconservative. 
As shown in Table 6.1, the incorrect wall force is predicted whenever the scheme is 
nonconservative. In the Method I case (see also Fig. 6.3.1), the phase of the wall force 
is in error by about 32.2®. The phase error in the Method II example is 6.1®. 

From these numerical results we conclude that both Method I and Method II 
linearizations will produce satisfactory results if and only if the linearizations are 
conservative. However, since the Method I linearization is predicated on the assump- 
tion that a constant time step is used throughout the computational domain, this 
precludes the use of Method I for most problems since it would be difficult and un- 
desirable to generate computational grids with constant cell areas throughout the 
computational domain. With the Method II analysis, we only require that AtfAA 
be constant for the scheme to be conservative. Therefore, for the remaining examples. 
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we will use a conservative Method II analysis. 

6.3.2 Unsteady Compressor and Fan Flows 

Having demonstrated the ability of the present method to model transonic channel 
flow, we next consider the unsteady flow in compressors and fans. 

Tenth Standard Configuration 

The first cascade considered is the Tenth Standard Configuration [9]. The airfoils of 
this cascade have a NACA 0006 thickness distribution slightly modified so that the 
trailing edge is wedged rather than blunt. The camber line is a circular arc with a 
maximum height of 5 percent of the chord. The flow conditions are such that there 
is a supersonic patch on the suction surface of the airfoil. The stagger angle, 0, is 
45* and the gap-to-chord ratio, G, is 1.0. The mean inflow angle, j3_^, is 58° and 
the inflow Mach number, M_^, is 0.8. Figure 6.3.2 shows the computed coefficient of 
pressure distribution along the airfoil surface calculated using the present nonlinear 
steady Euler code. The grid used for this calculation was a 193 x 49 node H-grid with 
a total of 193 nodes on the airfoil surface. Note in particular the transonic patch on 
the suction surface of the airfoil. The present steady Euler solver captures the shock 
over about five grid points. Also shown is the a nonlinear Euler solution provided by 
Huff based on a flux difference splitting algorithm [30,31]. 

With the steady solution now known, consider the case where the airfoils plunge 
with an interblade phase angle, a, of —90° and a reduced frequency, uJ (based on the 
upstream velocity and blade chord), of 1.287. Figure 6.3.2 shows the computed un- 
steady pressure distribution on the airfoil surface using Method II linearization. The 
impulsive shock load is clearly visible on the suction surface. Also shown for com- 
parison is the pressure distribution computed using Huff’s nonlinear time-marching 
algorithm. The agreement between the present linearized analysis and the nonlinear 
time-marching Euler analysis is excellent away from the shock. Shown in the table 
insert in Fig. 6.3.2 is the magniitude and phcise of the resulting unsteady lift. The 
magnitude of the unsteady lift calculated using the two different approaches agrees 
within about 2%; the phase differs by only about 3°. Note that the shock impulse 
predicted by the present unsteady linearized Euler analysis is somewhat narrower 
and taller than that predicted by the nonlinear code. The areas of the impulses, 
however, are very nearly equal. Furthermore, the unsteady load due to the impulse 
is of the same order of magnitude as the unsteady load due to the unsteady pressure 
distribution away from the shock. 

Because conservative Method II linearizations require that Af/AA be constant 
throughout the computational domain, the time step taken in a particular compu- 
tational cell may be considerably smaller than the maximum permitted for stable 
calculations. The result is that the convergence will be considerably slower than 
if the local maximum permissible time step had been taken everywhere (local time 
stepping). To overcome this problem, we use conservative time stepping in conjunc- 
tion with multiple grid acceleration. Figure 6.3.2 shows the convergence histories for 
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Figure 6.5: Coefficient of pressure distribution, Tenth Standard Configuration: 

= 0.8, G = 1.0, © = 45», = 58° 
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Figure 6.6: Real and imaginary unsteady surface pressure, Tenth Standard Configu- 
ration, plunging: u; = 1.287, = —90® 
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Iterations 

Figure 6.7: Convergence histories of unsteady solution for different methods. All 
cases use multigrid acceleration. 

three linearized unsteady flow calculations for the previous example: one using local 
time stepping with multigrid, one using conservative time stepping with multigrid, 
and one using conservative time stepping with over-relaxation plus multigrid. Note 
that over-relaxation reduces the computational time required by a factor of about 2.5 
compared to conservative time stepping without over-relaxtion. Finally, we should 
mention that a comparable nonlinear time marching algorithm would require about 
20 to 50 times the computational time required by the global time step calculations 
with over-relaxation and multigrid. 

High Speed Cascade 

The next case considered is a two-dimensional cascade of fan blades with a relative 
inlet Mach number, M.^o, of 1.2, stagger angle, 0, of 55®, and blade-to-blade gap, G, 
of 1.0. This case is presented to demonstrate the importeince of moving shocks on the 
aeroeleistic response of fau blades. Figure 6.3.2 shows the steady pressure contours. 
The solution weis computed on a 129 x 33 node grid with a total of 129 nodes on 
the airfoil surface. Figure 6.3.2 shows the computed isentropic Mach number on the 
airfoil’s surface. The pressure rise due to the passage shock can be clearly seen on 
both the suction and pressure surfaces. The shock is smeared over about four grid 
nodes. 
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Figure 6,8: Steady pressure contours, modified circular arc airfoil: = 1.2, G = 

l.O,0 = 55“,f2_^ = 60» 
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Figure 6.9: Isentropic Mach number distribution, modified circular arc airfoil: Af_, 
1.2,G = l.O,0 = 55®,n_oo = 6O° 
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Interblade Phase Angle, ct 

Figure 6.10: Imaginary part of moment coefficient for a range of interblade phase 
angles, modified circular arc airfoil, pitching about midchord, a; = 0.5. 


Next, we computed the unsteady aerodynamic response of the cascade for a range 
of interblade phase angles. The airfoils pitch about their midchords with a reduced 
frequency, U, of 0.5. For each interblade phase angle, the computed unsteady surface 
pressure was integrated to obtain the unsteady pitching moment. Shown in Fig 6.3.2 is 
the imaginary part of the unsteady pitching moment as a function of interblaxle phase 
angle. Positive imaginary pitching moments correspond to negative aerodynamic 
damping which will produce flutter for tuned cascades. Note that for er = 120°, the 
cascade is unstable. 

Shown in Fig. 6.3.2 is the unsteady pressure for the case where the airfoils vibrate 
in pitch with a reduced frequency, uJ, of 0.5 and interblade phase angle, <t, of 120°. 
Note that the unsteady aerodynamic load on the airfoil is dominated by the shock 
impulses. The impulse acting near the trailing edge provides a positive contribution 
to the imaginary part and hence is destabilizing. The impulse near the leading edge, 
on the other hand, is stabilizing. 

While these results demonstrate the importance of unsteady shock motion on the 
unsteady aerodynamic behavior of the fan, it should be noted that whenever strong 
in-passage shocks occur, viscous effects become important due to the large adverse 
pressure gradient at the shock. The actual unsteady aerodynamic behavior may be 
substantially different. 
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Figure 6.11: Real and imaginary unsteady surface pressure, modified circular arc 
airfoil, pitching about midchord: u = 0.5, <t = 120® 
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6.4 Summary 

In this chapter, we have presented a linearized Euler analysis of two-dimensional un- 
steady transonic flows in channels and cascades. Two different types of linearization 
were examined. Using Method I, the nonlinear Euler equations are first discretized 
using a conservative, time-accurate Lax-Wendroff scheme. The resulting nonlinear 
finite volume discretization is then linearized. Using Method II, the Euler equations 
are first linearized and then discretized using a Lax-Wendroff scheme. It was shown 
mathematically and by numerical experiment that the both Method I and Method 
II linearizations correctly predict the unsteady shock impulse in transonic flows if 
and only if the scheme is conservative; the order of linearization appears to be in- 
consequential. When either the Method I or Method II discretizations were made 
nonconservative by using a non-constant At/ A A, the shock impulse was found to 
be incorrectly predicted even though the methods are formally second-order accurate 
and consistent with the linearized Euler equations. 

Because a constant At/AA is required in the steady and unsteady flow calcu- 
lations to insure conservation, the time step taken at a computational cell may be 
significantly smaller than the maximum local permissible time step. This small time 
step in turn slows convergence of the scheme. To overcome this difficulty, an over- 
relaxation technique was proposed that dramatically improves the convergence rate 
of the linearized Euler analysis while leaving the method fully conservative. When 
coupled with Ni’s multiple grid acceleration technique, the present linearized Euler 
solver can compute unsteady transonic flows nearly two orders-of-magnitude faster 
than a comparable nonlinear time-accurate time-marching solver. 

A number of two- and three-dimensional unsteady transonic flows in cascades 
were computed using the linearized Euler analyses. Where possible, these results 
were compared to a nonlinear time-accurate time-marching scheme and found to be 
in excellent agreement. Furthermore, the unsteady shock load was found to be a 
significant contributor to the unsteady aerodynamic forces acting on the airfoil. 
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Conclusions And Future 
Considerations 

7.1 Conclusions 

The main goal of the present research was to develop an accurate and efficient com- 
putational analysis that will enable aeroelasticians to understand unsteady aerody- 
namic phenomena encountered in turbomachinery. Toward that end, a linearized 
Euler solver capable of accurately determining the unsteady loads in cascades oper- 
ating in strongly transonic flow regimes at moderate reduced frequencies has been 
developed. The technique selected for both the present steady and unsteady analyses 
is bzLsed on Ni’s variation of the Lax-Wendroff scheme. A pseudo time dependence is 
introduced into the governing equations making them hyperbolic in time. Therefore, 
the resulting discretized equations can be marched in time using traditional CFD 
techniques. Since a steady-state solution is desired in both the present steady and 
unsteady analyses, the calculations are not restricted to time-accurate time marching 
and acceleration techniques such eis local time stepping and multiple grid accelera- 
tion are used. The result is the method requires one to two orders of magnitude 
less computational time than traditional nonlinear time-accurate time-marching al- 
gorithms. In the course of development, the present method has incorporated three 
major improvements not previously encountered in linearized analyses. These are: 

(1) A continuously deforming computational grid is used which is capable of con- 
forming to both rigid body and flexible blade motions. The use of a deforming grid 
eliminates error producing extrapolation terms that would otherwise appear in the 
flow tangency (solid surface) boundary condition, substantially increasing the accu- 
racy of the method. 

(2) Numerically exact nonreflecting boundary conditions were developed that elim- 
inate all spurious reflections of outgoing pressure, entropy and vorticity waves in the 
far fleld. The new boundary condition formulation is quite general and can be applied 
to other flow models (such as the potential euad Navier-Stokes equations) and can be 
extended to handle fully three-dimensional flows. 

(3) Numerical computations demonstrated the ability of the present method to 
predict accurately the unsteady loads associated with transonic operating conditions 
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using a linearized shock capturing technique. Until recently, many researchers be- 
lieved transonic unsteady flows to be fundamentally nonlinear, violating traditional 
small disturbance assumptions. In particular, questions have arisen as to whether 
strong in-passage shocks could be modeled within the framework of a linearized anal- 
ysis. The present research has demonstrated the ability to predict accurately the 
unsteady loading in such cases provided a conservative linearization is used. 

7.2 Future Considerations 

There are three issues left to address before the present two-dimensional linearized 
Euler analysis is complete: 

(1) Quasi-three-dimensional effects, due to variations of streamtube height, need 
to be incorporated into the present analysis. Modeling of realistic compressor and 
turbine stages requires the capability of specifying changes in streamtube height. 
Although the present method still will not be capable of analyzing fully three- dimen- 
sional geometries, the additional physics being modeled will improve the usefulness 
of the present analysis as a design tool. 

(2) A more efficient algorithm needs to be developed to compute the eigenmodes of 

the linearized unsteady flow solver in the far field. The application of the numerically 
exact far-field boundary requires the solution of a sparse generalized complex eigen- 
value problem. The present analysis uses a standard eigenvalue solver (EISPACK) 
which does not take advantage of the spaxseness of the far-field matrices. Since the 
far-field eigenanalysis currently requires approximately 25% of the total computa- 
tional time (for a 65 x 17 node grid) significant gains in speed can be obtained if a 
more efficient eigensolver were used that takes advantage of sparseness and the known 
analytical eigenmode shapes. , 

(3) To further validate the capabilities of the present linearized method to predict 
unsteady loads in cascades, additional comparisons to experimental data are required. 
Although some comparisons have been made and the results are very encouraging, 
the extent to which the code can accurately predict unsteady loading still needs to 
be quantified. 

Finally, in future work viscous effects need to be modeled within a linearized 
framework. Viscous effects are essential if one is to model transonic flow problems 
with passage shocks since viscous forces help fix the location of the shock. Addition- 
ally, flow problems where separation occurs, as in the case of stall flutter, cannot be 
modeled within the framework of inviscid analyses. It should be possible to extend 
the linearized approach developed in this research to the Navier- Stokes equations and 
to inviscid/ viscous interaction models. 
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